Testing Lanchester’s laws in coral reef fish
David Cerny
1/28/2018
Data preparation
Read in the data spreadsheets:
sheet <- read.csv("/Users/David/Documents/College/2017–18/Winter 2018/FBQ_git/LiveFishInteraction_Feb9.csv", sep = ",", strip.white = TRUE)
scored <- read.csv("/Users/David/Documents/College/2017–18/Winter 2018/FBQ_git/ScoredFishInteraction_Feb19.csv", sep = ",", strip.white = TRUE)
# Combining the two is easier when dealing with the relevant columns alone. These are:
# Species.A, Species.B, Group.A, Group.B, and Winner.
ABwin_live <- sheet[,c(4,6,8,9,10)]
ABwin_scored <- scored[,c(10,12,13,15,16)]
both <- rbind(ABwin_live, ABwin_scored)
# Create new versions of the three data frames that include one-on-one interactions only:
sheet_oneonone <- sheet[sheet$Group.A == 1 & sheet$Group.B == 1,]
scored_oneonone <- scored[scored$Group.A == 1 & scored$Group.B == 1,]
both_oneonone <- both[both$Group.A == 1 & both$Group.B == 1,]Clean-up
Step 1. Unidentified species are marked with an asterisk. Exclude them from the matrix based on this criterion:
drop_ufos <- function(df) {
# Dummy vector for storing the numbers of the rows with unidentifiable species:
tobedropped <- vector()
for (i in 1:nrow(df)) {
if (grepl("*", df$Species.A[i], fixed = TRUE) == TRUE |
grepl("*", df$Species.B[i], fixed = TRUE) == TRUE |
grepl("*", df$Winner[i], fixed = TRUE) == TRUE) {
tobedropped <- c(tobedropped, i)
}
}
return(df[-c(tobedropped),])
}
sheet2 <- drop_ufos(sheet)
scored2 <- drop_ufos(scored)
both2 <- drop_ufos(both)
sheet_ooo2 <- drop_ufos(sheet_oneonone)
scored_ooo2 <- drop_ufos(scored_oneonone)
both_ooo2 <- drop_ufos(both_oneonone)Step 2. Remove intraspecific interactions (a few of these were scored for live observations in the beginning, before it was decided to selectively focus on interspecific interactions only):
drop_intra <- function(df) {
tobedropped <- vector()
for (i in 1:nrow(df)) {
if (identical(as.character(df$Species.A[i]), as.character(df$Species.B[i])) == TRUE) { tobedropped <- c(tobedropped, i)
}
}
if (length(tobedropped) > 0) {
return(df[-c(tobedropped),])
} else {
return(df)
}
}
sheet3 <- drop_intra(sheet2)
scored3 <- drop_intra(scored2)
both3 <- drop_intra(both2)
sheet_ooo3 <- drop_intra(sheet_ooo2)
scored_ooo3 <- drop_intra(scored_ooo2)
both_ooo3 <- drop_intra(both_ooo2)Step 3. Drop “Striated Surgeonfish” from the matrices, as this is likely a composite category comprising both striated and brown surgeonfish. These species may have different dominance ranks, but are almost indistinguishable in appearance.
no_surgeon <- function(df) {
return(df[df$Species.A != "Striated Surgeonfish" & df$Species.B != "Striated Surgeonfish",])
}
sheet4 <- no_surgeon(sheet3)
scored4 <- no_surgeon(scored3)
both4 <- no_surgeon(both3)
sheet_ooo4 <- no_surgeon(sheet_ooo3)
scored_ooo4 <- no_surgeon(scored_ooo3)
both_ooo4 <- no_surgeon(both_ooo3)Generating a dominance matrix
Extract the levels of Species.A and Species.B and use them as the rows and columns of the dominance matrix:
# The union of the species listed in the "Species A" and "Species B" columns will be used
# as the rows and columns of a new matrix filled with zeros:
emptymatrix <- function(df) {
species <- union(levels(factor(df$Species.A)), levels(factor(df$Species.B)))
rownum <- length(species)
newmatrix <- matrix(0, nrow = rownum, ncol = rownum)
row.names(newmatrix) <- species
colnames(newmatrix) <- species
return(newmatrix)
}
dommatrix <- emptymatrix(sheet4)
scoredmatrix <- emptymatrix(scored4)
bothmatrix <- emptymatrix(both4)
domooo <- emptymatrix(sheet_ooo4)
scoredooo <- emptymatrix(scored_ooo4)
bothooo <- emptymatrix(both_ooo4)
# If species A won over species B (winner = A), add 1 to the cell in column A and row B.
# If species B won over species A (winner = B), add 1 to the cell in column B and row A.
# An increment function will be used to increase the value of each cell accordingly:
inc <- function(x) # credit: Grega Kešpret, Stack Overflow
{
eval.parent(substitute(x <- x + 1))
}
# Now apply it to populate the matrix:
populate <- function(sheet, df) {
for (i in 1:nrow(sheet)) {
spec_a <- as.character(sheet$Species.A[i])
spec_b <- as.character(sheet$Species.B[i])
winner <- as.character(sheet$Winner[i])
if (identical(spec_a, winner) == TRUE) {
inc(df[as.character(spec_b), as.character(spec_a)])
}
if (identical(spec_b, winner) == TRUE) { # This should be equivalent to "else"
inc(df[as.character(spec_a), as.character(spec_b)])
}
}
diag(df) <- NA
return(df)
}
dommatrix <- populate(sheet4, dommatrix)
scoredmatrix <- populate(scored4, scoredmatrix)
bothmatrix <- populate(both4, bothmatrix)
domooo <- populate(sheet_ooo4, domooo)
scoredooo <- populate(scored_ooo4, scoredooo)
bothooo <- populate(both_ooo4, bothooo)Convert from matrices to data frames for easier manipulation:
convert_to_df <- function(df) {
frame <- as.data.frame(df, stringsAsFactors = FALSE)
for (i in 1:ncol(df)) {
frame[,i] <- as.numeric(frame[,i])
}
return(frame)
}
domall <- convert_to_df(dommatrix)
scoredall <- convert_to_df(scoredmatrix)
bothall <- convert_to_df(bothmatrix)
liveoooall <- convert_to_df(domooo)
scoredoooall <- convert_to_df(scoredooo)
bothoooall <- convert_to_df(bothooo)Exclude cornetfish, the only predator in the dataset. (Note: the fact that this is done here, at the matrix formatting stage rather than at the data clean-up stage, simply reflects the point at which the decision was taken during the actual analysis. Excluding a taxon from the matrix using the method below is equivalent to excluding it from the data sheet.)
# Exclude a taxon from a matrix:
# Original unvectorized function used in the pre-revisions part of the file is commented
# out below:
#
# exclude <- function(df, taxon) {
# colrow <- c(which(row.names(df) == taxon), which(colnames(df) == taxon))
# excluded <- df[-colrow[1], -colrow[2]]
# return(excluded)
# }
exclude <- function(df, taxon) {
retained <- df[!rownames(df) %in% taxon, !colnames(df) %in% taxon]
return(retained)
}
# Note: since there are no interactions involving cornetfish in the video-only data,
# we only need to exclude it from the live-only and combined datasets.
nocornet_domall <- exclude(domall, "Cornetfish")
nocornet_bothall <- exclude(bothall, "Cornetfish")
nocornet_live <- exclude(liveoooall, "Cornetfish")
nocornet_both <- exclude(bothoooall, "Cornetfish")Matrix completeness:
compl <- function(df) {
counter <- 0
for (i in 1:nrow(df)) {
for (j in 1:ncol(df)) {
if (!is.na(df[i,j]) & df[i,j] == 0 & df[j,i] == 0) {
counter <- counter + 1
}
}
}
return(100*(1 - (counter/(nrow(df)*(nrow(df) - 1)))))
}
paste("The completeness of the live-only all-species matrix is ", compl(nocornet_domall), "%", sep = ""); paste("The completeness of the video-only all-species matrix is ", compl(scoredall), "%", sep = ""); paste("The completeness of the combined all-species matrix is ", compl(nocornet_bothall), "%", sep = ""); paste("The completeness of the live-only one-on-one all-species matrix is ", compl(nocornet_live), "%", sep = ""); paste("The completeness of the video-only one-on-one all-species matrix is ", compl(scoredoooall), "%", sep = ""); paste("The completeness of the combined one-on-one all-species matrix is ", compl(nocornet_both), "%", sep = "")## [1] "The completeness of the live-only all-species matrix is 38.961038961039%"## [1] "The completeness of the video-only all-species matrix is 27.2058823529412%"## [1] "The completeness of the combined all-species matrix is 37.9446640316206%"## [1] "The completeness of the live-only one-on-one all-species matrix is 33.7662337662338%"## [1] "The completeness of the video-only one-on-one all-species matrix is 25%"## [1] "The completeness of the combined one-on-one all-species matrix is 33.596837944664%"Landau’s h-index:
library(gdata)
library(compete)
# To determine linearity, we need to make our dominance matrix symmetric. Below, we do
# this by subtracting the lower triangle from the upper triangle, and flipping the results
# around the diagonal:
symmatrix <- function(df) {
# Going through the upper triangle column-wise corresponds to going through the lower
# triangle row-wise:
triangle <- as.numeric(upperTriangle(df)) - as.numeric(lowerTriangle(df), byrow = TRUE)
# Replace the upper triangle in the default column-wise order:
symmat <- matrix(NA, nrow = nrow(df), ncol = ncol(df))
upperTriangle(symmat) <- triangle
# Replace the lower triangle in the row-wise order:
lowerTriangle(symmat, byrow = TRUE) <- triangle
return(symmat)
}
# Note that the p-values are simulation-based, and therefore stochastic: they will differ
# slightly between runs
devries(nocornet_domall, Nperms = 10000, history = FALSE, plot = TRUE)## h-modified = 0.3196588
## p-value= 0.0015
devries(scoredall, Nperms = 10000, history = FALSE, plot = TRUE)## h-modified = 0.2982113
## p-value= 0.0528
devries(nocornet_bothall, Nperms = 10000, history = FALSE, plot = TRUE)## h-modified = 0.2962757
## p-value= 0.002
devries(nocornet_live, Nperms = 10000, history = FALSE, plot = TRUE)## h-modified = 0.2779799
## p-value= 0.0097
devries(scoredoooall, Nperms = 10000, history = FALSE, plot = TRUE)## h-modified = 0.2972647
## p-value= 0.0559
devries(nocornet_both, Nperms = 10000, history = FALSE, plot = TRUE)## h-modified = 0.2678731
## p-value= 0.0066
Uncomment the code below (and change the paths as needed) to print the matrices:
# write.table(nocornet_domall, "/Users/David/Documents/College/2017–18/Winter 2018/FBQ_git/live_only_matrix.txt", quote = FALSE, sep = "\t")
# write.table(scoredall, "/Users/David/Documents/College/2017–18/Winter 2018/FBQ_git/video_only_matrix.txt", quote = FALSE, sep = "\t")
# write.table(nocornet_bothall, "/Users/David/Documents/College/2017–18/Winter 2018/FBQ_git/both_datasets_matrix.txt", quote = FALSE, sep = "\t")
# write.table(nocornet_live, "/Users/David/Documents/College/2017–18/Winter 2018/FBQ_git/live_only_one-on-one_matrix.txt", quote = FALSE, sep = "\t")
# write.table(scoredoooall, "/Users/David/Documents/College/2017–18/Winter 2018/FBQ_git/video_only_one-on-one_matrix.txt", quote = FALSE, sep = "\t")
# write.table(nocornet_both, "/Users/David/Documents/College/2017–18/Winter 2018/FBQ_git/both_datasets_one-on-one_matrix.txt", quote = FALSE, sep = "\t")CBI
The code below was used to compute the Clutton-Brock et al. dominance index (Clutton-Brock et al. 1979, 1982).
# Remember:
# If A won over B, add 1 to [B,A]
# If A won over C, add 1 ro [C,A]
# -> So the A column represents the species that A dominates (A wins)
# If A lost to B, add 1 to [A,B]
# If A lost to C, add 1 to [A,C]
# -> So the A row represents the species that dominate A (A losses)
# (B + b + 1)/(L + l + 1)
cbi <- function(df, species) {
# number of entries in column 'species' that are neither 0 nor NA:
B <- sum(df[,species] != 0, na.rm = TRUE)
# species corresponding to those entries:
dominates <- rownames(df)[which(df[,species] != 0)]
# number of entries in the columns corresponding to those species that are neither 0 nor NA:
if (length(dominates) > 0) {
b <- sum(subset(df, select = dominates) != 0, na.rm = TRUE)
} else {
b <- 0
}
# number of entries in row B that are neither 0 nor NA:
L <- sum(df[species,], na.rm = TRUE)
# species corresponding to those entries:
dominated <- colnames(df)[which(df[species,] != 0)]
# number of entries in the rows corresponding to those species that are neither 0 nor NA:
if (length(dominated) > 0) {
l <- sum(subset(df, rownames(df) %in% dominated) != 0, na.rm = TRUE)
} else {
l <- 0
}
return((B + b + 1)/(L + l + 1))
}
# Apply it to the full matrix:
cbi_hierarchy <- function(df) {
cbi_matrix <- data.frame(matrix(NA, nrow = nrow(df), ncol = 2))
colnames(cbi_matrix) <- c("Species", "CBI")
for (i in 1:nrow(df)) {
cbi_matrix$Species[i] <- rownames(df)[i]
cbi_matrix$CBI[i] <- cbi(df, rownames(df)[i])
}
cbi_matrix <- cbi_matrix[order(cbi_matrix$CBI, decreasing = TRUE),]
return(cbi_matrix)
}
cbi_ncdomall <- cbi_hierarchy(nocornet_domall)
cbi_scoredall <- cbi_hierarchy(scoredall)
cbi_ncbothall <- cbi_hierarchy(nocornet_bothall)
cbi_ncooolive <- cbi_hierarchy(nocornet_live)
cbi_ooovideo <- cbi_hierarchy(scoredoooall)
cbi_ncoooboth <- cbi_hierarchy(nocornet_both)Frequency of winning
win_freq <- function(df) {
j <- 1
species <- union(levels(df$Species.A), levels(df$Species.B))
freq_table <- data.frame(matrix(NA, nrow = length(species), ncol = 2))
colnames(freq_table) <- c("Species", "WinFreq")
for(i in species) {
if (grepl("*", i, fixed = TRUE) == TRUE) {
next
} else {
freq_table[j,1] <- i
involved <- sum(df$Species.A == i) + sum(df$Species.B == i)
won <- sum(df$Winner == i)
freq_table[j,2] <- won/involved
j <- j + 1
}
}
return(freq_table[complete.cases(freq_table),])
}Sizes
Compare the sizes estimated from live observations with those reported by Randall (2005):
# Maximum attained sizes according to Randall (2005):
sizes <- read.csv("/Users/David/Documents/College/2017–18/Winter 2018/FBQ_git/FishSize.csv", sep = ";", strip.white = TRUE, header = TRUE)
colnames(sizes)[1] <- "Species"
booksizes <- sizes$Attain[c(1:3,5:8,10:24)]
# Distributions of sizes observed in the field:
fishsize <- read.csv("/Users/David/Documents/College/2017–18/Winter 2018/FBQ_git/FieldFishSize.csv", sep = ",", strip.white = TRUE)
# Calculate the mean, median, and maximum of the observed size distributions:
meansizes <- colMeans(fishsize[,c(2:4,6:9,11:25)], na.rm = TRUE)
colMedians <- function(matrix) {
return(apply(matrix, 2, median, na.rm = TRUE))
}
colMax <- function(matrix) {
return(apply(matrix, 2, max, na.rm = TRUE))
}
mediansizes <- colMedians(fishsize[,c(2:4,6:9,11:25)])
maxsizes <- colMax(fishsize[,c(2:4,6:9,11:25)])
# Do book sizes correlate with live observation sizes?
# a <- cor.test(meansizes, mediansizes, method = "spearman", alternative = "two.sided")
# paste("correlation: ", signif(a$estimate[[1]],5), ", p-value: ", a$p.value, sep = "")
# b <- cor.test(meansizes, maxsizes, method = "spearman", alternative = "two.sided")
# paste("correlation: ", signif(b$estimate[[1]],5), ", p-value: ", b$p.value, sep = "")
# c <- cor.test(mediansizes, maxsizes, method = "spearman", alternative = "two.sided")
# paste("correlation: ", signif(c$estimate[[1]],5), ", p-value: ", c$p.value, sep = "")
d <- cor.test(meansizes, mediansizes, method = "pearson", alternative = "two.sided")
paste("correlation: ", signif(d$estimate[[1]],5), ", p-value: ", d$p.value, sep = "")## [1] "correlation: 0.9969, p-value: 1.4440712484574e-23"e <- cor.test(meansizes, maxsizes, method = "pearson", alternative = "two.sided")
paste("correlation: ", signif(e$estimate[[1]],5), ", p-value: ", e$p.value, sep = "")## [1] "correlation: 0.94797, p-value: 2.11724427906556e-11"f <- cor.test(mediansizes, maxsizes, method = "pearson", alternative = "two.sided")
paste("correlation: ", signif(f$estimate[[1]],5), ", p-value: ", f$p.value, sep = "")## [1] "correlation: 0.94834, p-value: 1.973173161907e-11"# Non-unique ranks, gives an error:
# g <- cor.test(booksizes, mediansizes, method = "spearman", alternative = "two.sided")
# paste("correlation: ", signif(g$estimate[[1]],5), ", p-value: ", g$p.value, sep = "")
# h <- cor.test(booksizes, meansizes, method = "spearman", alternative = "two.sided")
# paste("correlation: ", signif(h$estimate[[1]],5), ", p-value: ", h$p.value, sep = "")
# i <- cor.test(booksizes, maxsizes, method = "spearman", alternative = "two.sided")
# paste("correlation: ", signif(i$estimate[[1]],5), ", p-value: ", i$p.value, sep = "")
j <- cor.test(booksizes, mediansizes, method = "pearson", alternative = "two.sided")
paste("correlation: ", signif(j$estimate[[1]],5), ", p-value: ", j$p.value, sep = "")## [1] "correlation: 0.83931, p-value: 1.04831251950256e-06"k <- cor.test(booksizes, meansizes, method = "pearson", alternative = "two.sided")
paste("correlation: ", signif(k$estimate[[1]],5), ", p-value: ", k$p.value, sep = "")## [1] "correlation: 0.85355, p-value: 4.41082818979503e-07"l <- cor.test(booksizes, maxsizes, method = "pearson", alternative = "two.sided")
paste("correlation: ", signif(l$estimate[[1]],5), ", p-value: ", l$p.value, sep = "")## [1] "correlation: 0.83313, p-value: 1.48832376900543e-06"# Reassign the mediansizes vector so that it includes Juvenile Parrotfish:
mediansizes <- colMedians(fishsize[,c(2:4,6:25)])
# Create a new data frame that associates each species with the median of its observed
# size distribution, excluding Cornetfish and "Striated Surgeonfish":
obs_sizes <- data.frame(matrix(ncol = 2, nrow = length(mediansizes)))
colnames(obs_sizes) <- c("Species", "Size")
obs_sizes[,1] <- as.character(sizes$Species[c(1:3,5:24)])
obs_sizes[,2] <- as.numeric(mediansizes)Include mass (scalled isometrically and allometrically) as an additional variable:
mass <- read.csv("/Users/David/Documents/College/2017–18/Winter 2018/FBQ_git/body_to_mass_data.csv", sep = ",", strip.white = TRUE)
colnames(mass)[1] <- "Species"
lengthmass <- merge(obs_sizes, mass, by = "Species")
lengthmass[,5] <- lengthmass[,3]*(lengthmass[,2]^lengthmass[,4])
lengthmass[,6] <- lengthmass[,5]^(2/3)
colnames(lengthmass)[c(2,5,6)] <- c("Length", "Mass", "Mass2_3")Test of mass vs. individual fighting ability
Test the relationship between the two measures of dominance (CBI and winning frequency) and body mass or body mass raised to 2/3 by creating a single data frame that contains all the four variables. Note: from this point on, we focus exclusively on the combined (live plus video) dataset of 23 species (no cornetfish) and one-on-one-interactions.
both_test <- merge(cbi_ncoooboth, win_freq(both_ooo4[both_ooo4$Species.A != "Cornetfish" & both_ooo4$Species.B != "Cornetfish",]), by = "Species")
both_test <- merge(both_test, lengthmass, by = "Species")
# Number of interactions:
nrow(both_ooo4[both_ooo4$Species.A != "Cornetfish" & both_ooo4$Species.B != "Cornetfish",])## [1] 401cbi_mass <- lm(both_test$CBI ~ both_test$Mass)
summary(cbi_mass)$coefficients[2, 4]## [1] 0.4754628summary(cbi_mass)$r.squared## [1] 0.02452714winfreq_mass <- lm(both_test$WinFreq ~ both_test$Mass)
summary(winfreq_mass)$coefficients[2, 4]## [1] 0.1346416summary(winfreq_mass)$r.squared## [1] 0.1033787cbi_mass2_3 <- lm(both_test$CBI ~ both_test$Mass2_3)
summary(cbi_mass2_3)$coefficients[2, 4]## [1] 0.2865715summary(cbi_mass2_3)$r.squared## [1] 0.05387003winfreq_mass2_3 <- lm(both_test$WinFreq ~ both_test$Mass2_3)
summary(winfreq_mass2_3)$coefficients[2, 4]## [1] 0.1029628summary(winfreq_mass2_3)$r.squared## [1] 0.1215855# How far is spotted porcupinefish from the median in terms of the interquartile range?
# (Note: this is a simple measure of the extremeness of an outlier)
(both_test$Mass[both_test$Species == "Spotted Porcupinefish"] - median(both_test$Mass))/IQR(both_test$Mass)## [1] 9.456649(both_test$Mass2_3[both_test$Species == "Spotted Porcupinefish"] - median(both_test$Mass2_3))/IQR(both_test$Mass2_3)## [1] 4.972994# Exclude spotted porcupinefish:
nocornet_nosp_both <- exclude(nocornet_both, "Spotted Porcupinefish")
cbi_ncoooboth_nosp <- cbi_hierarchy(nocornet_nosp_both)
both_test_nosp <- merge(cbi_ncoooboth_nosp, win_freq(both_ooo4[!both_ooo4$Species.A %in% c("Cornetfish", "Spotted Porcupinefish") & !both_ooo4$Species.B %in% c("Cornetfish", "Spotted Porcupinefish"),]), by = "Species")
both_test_nosp <- merge(both_test_nosp, lengthmass, by = "Species")
# Number of data points:
nrow(both_ooo4[!both_ooo4$Species.A %in% c("Cornetfish", "Spotted Porcupinefish") & !both_ooo4$Species.B %in% c("Cornetfish", "Spotted Porcupinefish"),])## [1] 398cbi_mass_nosp <- lm(both_test_nosp$CBI ~ both_test_nosp$Mass)
summary(cbi_mass_nosp)$coefficients[2, 4]## [1] 0.0301074summary(cbi_mass_nosp)$r.squared## [1] 0.214146winfreq_mass_nosp <- lm(both_test_nosp$WinFreq ~ both_test_nosp$Mass)
summary(winfreq_mass_nosp)$coefficients[2, 4]## [1] 0.2470032summary(winfreq_mass_nosp)$r.squared## [1] 0.06638793cbi_mass2_3_nosp <- lm(both_test_nosp$CBI ~ both_test_nosp$Mass2_3)
summary(cbi_mass2_3_nosp)$coefficients[2, 4]## [1] 0.03673997summary(cbi_mass2_3_nosp)$r.squared## [1] 0.2003125winfreq_mass2_3_nosp <- lm(both_test_nosp$WinFreq ~ both_test_nosp$Mass2_3)
summary(winfreq_mass2_3_nosp)$coefficients[2, 4]## [1] 0.2267295summary(winfreq_mass2_3_nosp)$r.squared## [1] 0.07216368This is the code used to generate Figure 1 of the manuscript (NOTE: This figure summarizes the results of the regression analyses uncorrected for phylogenetic relatedness, and is now given as Figure S3 in the Supplementary Information.)
# png("/Users/David/Documents/College/2017–18/Winter 2018/FBQ_git/Figure_1.png", width = 3600, height = 3600, pointsize = 96)
png("/Users/David/Documents/College/2017–18/Winter 2018/FBQ_git/Supp_Figure_3.png", width = 3600, height = 3600, pointsize = 96)
# setEPS()
# postscript("/Users/David/Documents/College/2017–18/Winter 2018/FBQ_git/Figure_1.eps")
spec_dec <- function(x, k) trimws(format(round(x, k), nsmall = k)) # credit: Jeromy Anglim, Stack Overflow
layout(matrix(c(seq(1,4,1)), nrow=2, ncol=2, byrow=TRUE))
par(mar=c(5, 4, 2, 2))
plot(both_test$Mass, both_test$CBI, xlab = "Body mass (g)", ylab = "CBI", main = "Isometric scaling", pch = 16, col = ifelse(both_test$Mass == both_test$Mass[both_test$Species == "Spotted Porcupinefish"], "gray70", "black"))
title("(a)", line = -1, adj = 0, outer = TRUE)
title("(b)", line = -1, outer = TRUE)
title("(c)", line = -20, adj = 0, outer = TRUE)
title("(d)", line = -20, outer = TRUE)
abline(lm(both_test$CBI ~ both_test$Mass), lwd = 5)
abline(lm(both_test_nosp$CBI ~ both_test_nosp$Mass), lwd = 5, col = "gray70")
rsquared <- summary(cbi_mass)$r.squared
pvalue <- summary(cbi_mass)$coefficients[2,4]
rsquared_nosp <- summary(cbi_mass_nosp)$r.squared
pvalue_nosp <- summary(cbi_mass_nosp)$coefficients[2,4]
text(800, 10.4, cex = 0.8, bquote(italic(R)^2 == .(spec_dec(rsquared, 3))), pos = 4)
text(800, 9.1, cex = 0.8, bquote(italic(p) == .(spec_dec(pvalue, 3))), pos = 4)
text(800, 7.3, cex = 0.8, bquote(italic(R)^2 == .(spec_dec(rsquared_nosp, 3))), col = "gray70", pos = 4)
text(800, 6, cex = 0.8, bquote(italic(p) == .(spec_dec(pvalue_nosp, 3))), col = "gray70", pos = 4)
par(mar=c(5, 4, 2, 2))
plot(both_test$Mass2_3, both_test$CBI, xlab = "Body mass raised to 2/3", ylab = "", main = "Allometric scaling", pch = 16, col = ifelse(both_test$Mass == both_test$Mass[both_test$Species == "Spotted Porcupinefish"], "gray70", "black"))
abline(lm(both_test$CBI ~ both_test$Mass2_3), lwd = 5)
abline(lm(both_test_nosp$CBI ~ both_test_nosp$Mass2_3), lwd = 5, col = "gray70")
rsquared <- summary(cbi_mass2_3)$r.squared
pvalue <- summary(cbi_mass2_3)$coefficients[2,4]
rsquared_nosp <- summary(cbi_mass2_3_nosp)$r.squared
pvalue_nosp <- summary(cbi_mass2_3_nosp)$coefficients[2,4]
text(15, 10.4, cex = 0.8, bquote(italic(R)^2 == .(spec_dec(rsquared, 3))), pos = 4)
text(15, 9.1, cex = 0.8, bquote(italic(p) == .(spec_dec(pvalue, 3))), pos = 4)
text(15, 7.3, cex = 0.8, bquote(italic(R)^2 == .(spec_dec(rsquared_nosp, 3))), col = "gray70", pos = 4)
text(15, 6, cex = 0.8, bquote(italic(p) == .(spec_dec(pvalue_nosp, 3))), col = "gray70", pos = 4)
par(mar=c(5, 4, 2, 2))
plot(both_test$Mass, both_test$WinFreq, xlab = "Body mass (g)", ylab = "Frequency of winning", pch = 16, col = ifelse(both_test$Mass == both_test$Mass[both_test$Species == "Spotted Porcupinefish"], "gray70", "black"))
abline(lm(both_test$WinFreq ~ both_test$Mass), lwd = 5)
abline(lm(both_test_nosp$WinFreq ~ both_test_nosp$Mass), lwd = 5, col = "gray70")
rsquared <- summary(winfreq_mass)$r.squared
pvalue <- summary(winfreq_mass)$coefficients[2,4]
rsquared_nosp <- summary(winfreq_mass_nosp)$r.squared
pvalue_nosp <- summary(winfreq_mass_nosp)$coefficients[2,4]
text(800, 0.4, cex = 0.8, bquote(italic(R)^2 == .(spec_dec(rsquared, 3))), pos = 4)
text(800, 0.4 - 1.3/16, cex = 0.8, bquote(italic(p) == .(spec_dec(pvalue, 3))), pos = 4)
text(800, 0.4 - (1.3+1.8)/16, cex = 0.8, bquote(italic(R)^2 == .(spec_dec(rsquared_nosp, 3))), col = "gray70", pos = 4)
text(800, 0.4 - (2*1.3+1.8)/16, cex = 0.8, bquote(italic(p) == .(spec_dec(pvalue_nosp, 3))), col = "gray70", pos = 4)
par(mar=c(5, 4, 2, 2))
plot(both_test$Mass2_3, both_test$WinFreq, xlab = "Body mass raised to 2/3", ylab = "", pch = 16, col = ifelse(both_test$Mass == both_test$Mass[both_test$Species == "Spotted Porcupinefish"], "gray70", "black"))
abline(lm(both_test$WinFreq ~ both_test$Mass2_3), lwd = 5)
abline(lm(both_test_nosp$WinFreq ~ both_test_nosp$Mass2_3), lwd = 5, col = "gray70")
rsquared <- summary(winfreq_mass2_3)$r.squared
pvalue <- summary(winfreq_mass2_3)$coefficients[2,4]
rsquared_nosp <- summary(winfreq_mass2_3_nosp)$r.squared
pvalue_nosp <- summary(winfreq_mass2_3_nosp)$coefficients[2,4]
text(60, 0.4, cex = 0.8, bquote(italic(R)^2 == .(spec_dec(rsquared, 3))), pos = 4)
text(60, 0.4 - 1.3/16, cex = 0.8, bquote(italic(p) == .(spec_dec(pvalue, 3))), pos = 4)
text(60, 0.4 - (1.3+1.8)/16, cex = 0.8, bquote(italic(R)^2 == .(spec_dec(rsquared_nosp, 3))), col = "gray70", pos = 4)
text(60, 0.4 - (2*1.3+1.8)/16, cex = 0.8, bquote(italic(p) == .(spec_dec(pvalue_nosp, 3))), col = "gray70", pos = 4)
dev.off()Testing Lanchester’s Laws
For each interaction, multiply the body mass of either species (scaled isometrically or allometrically) by its group size or group size squared. See which species maximizes the product, and whether it matches the observed winner. Run a binomial test on the proportion of successful predictions.
both_group <- both4[both4$Group.A != both4$Group.B,]
# Clean-up: drop rows containing NAs:
both_group <- both_group[complete.cases(both_group),]
# Exclude cornetfish:
both_group <- both_group[both_group$Species.A != "Cornetfish" &
both_group$Species.B != "Cornetfish",]
# Number of data points:
nrow(both_group)## [1] 83# Get the body sizes for both competitors:
for(i in 1:nrow(both_group)) {
species_a <- as.character(both_group$Species.A[i])
species_b <- as.character(both_group$Species.B[i])
both_group[i,6] <- lengthmass$Mass[lengthmass$Species == species_a]
both_group[i,7] <- lengthmass$Mass[lengthmass$Species == species_b]
}
colnames(both_group)[6:7] <- c("BodySize.A", "BodySize.B")
# Convert group sizes to the right object type (from 'character' to 'double'):
both_group$Group.A <- as.numeric(both_group$Group.A)
both_group$Group.B <- as.numeric(both_group$Group.B)
# Make predictions based on the linear law:
for(i in 1:nrow(both_group)) {
if (both_group$Group.A[i]*both_group$BodySize.A[i] > both_group$Group.B[i]*both_group$BodySize.B[i]) {
both_group[i,8] <- both_group$Species.A[i]
} else {
both_group[i,8] <- both_group$Species.B[i]
}
}
# Make predictions based on the square law:
for(i in 1:nrow(both_group)) {
if ((both_group$Group.A[i]^2)*both_group$BodySize.A[i] > (both_group$Group.B[i]^2)*both_group$BodySize.B[i]) {
both_group[i,9] <- both_group$Species.A[i]
} else {
both_group[i,9] <- both_group$Species.B[i]
}
}
colnames(both_group)[8:9] <- c("Linear.Law.Prediction", "Square.Law.Prediction")
# Test if the linear law works:
for(i in 1:nrow(both_group)) {
if (identical(as.character(both_group$Winner[i]), as.character(both_group$Linear.Law.Prediction[i])) == TRUE) {
both_group[i,10] <- 1
} else {
both_group[i,10] <- 0
}
}
# Test if the square law works:
for(i in 1:nrow(both_group)) {
if (identical(as.character(both_group$Winner[i]), as.character(both_group$Square.Law.Prediction[i])) == TRUE) {
both_group[i,11] <- 1
} else {
both_group[i,11] <- 0
}
}
colnames(both_group)[10:11] <- c("Linear.True", "Square.True")
# Proportion of successful predictions:
sum(both_group$Linear.True)/nrow(both_group)## [1] 0.6626506sum(both_group$Square.True)/nrow(both_group)## [1] 0.4939759# Test significance:
dbinom(sum(both_group$Linear.True), size = nrow(both_group), prob = 1/2) ## [1] 0.001053886dbinom(sum(both_group$Square.True), size = nrow(both_group), prob = 1/2)## [1] 0.08679764# When do the two laws give different predictions?
dif <- both_group[both_group$Linear.Law.Prediction != both_group$Square.Law.Prediction,]
dbinom(sum(dif$Linear.True), size = nrow(dif), prob = 1/2) ## [1] 0.001087189# Now repeat all of this for body mass raised to two thirds:
both_group2 <- both4[both4$Group.A != both4$Group.B,]
both_group2 <- both_group2[complete.cases(both_group2),]
both_group2 <- both_group2[both_group2$Species.A != "Cornetfish" &
both_group2$Species.B != "Cornetfish",]
for(i in 1:nrow(both_group2)) {
species_a <- as.character(both_group2$Species.A[i])
species_b <- as.character(both_group2$Species.B[i])
both_group2[i,6] <- lengthmass$Mass2_3[lengthmass$Species == species_a]
both_group2[i,7] <- lengthmass$Mass2_3[lengthmass$Species == species_b]
}
colnames(both_group2)[6:7] <- c("BodySize.A", "BodySize.B")
both_group2$Group.A <- as.numeric(both_group2$Group.A)
both_group2$Group.B <- as.numeric(both_group2$Group.B)
for(i in 1:nrow(both_group2)) {
if (both_group2$Group.A[i]*both_group2$BodySize.A[i] > both_group2$Group.B[i]*both_group2$BodySize.B[i]) {
both_group2[i,8] <- both_group2$Species.A[i]
} else {
both_group2[i,8] <- both_group2$Species.B[i]
}
}
for(i in 1:nrow(both_group2)) {
if ((both_group2$Group.A[i]^2)*both_group2$BodySize.A[i] > (both_group2$Group.B[i]^2)*both_group2$BodySize.B[i]) {
both_group2[i,9] <- both_group2$Species.A[i]
} else {
both_group2[i,9] <- both_group2$Species.B[i]
}
}
colnames(both_group2)[8:9] <- c("Linear.Law.Prediction", "Square.Law.Prediction")
for(i in 1:nrow(both_group2)) {
if (identical(as.character(both_group2$Winner[i]), as.character(both_group2$Linear.Law.Prediction[i])) == TRUE) {
both_group2[i,10] <- 1
} else {
both_group2[i,10] <- 0
}
}
for(i in 1:nrow(both_group2)) {
if (identical(as.character(both_group2$Winner[i]), as.character(both_group2$Square.Law.Prediction[i])) == TRUE) {
both_group2[i,11] <- 1
} else {
both_group2[i,11] <- 0
}
}
colnames(both_group2)[10:11] <- c("Linear.True", "Square.True")
sum(both_group2$Linear.True)/nrow(both_group2)## [1] 0.5662651sum(both_group2$Square.True)/nrow(both_group2)## [1] 0.3614458dbinom(sum(both_group2$Linear.True), size = nrow(both_group2), prob = 1/2) ## [1] 0.04240454dbinom(sum(both_group2$Square.True), size = nrow(both_group2), prob = 1/2)## [1] 0.003597748dif2 <- both_group2[both_group2$Linear.Law.Prediction != both_group2$Square.Law.Prediction,]
dbinom(sum(dif2$Linear.True), size = nrow(dif2), prob = 1/2) ## [1] 0.0001001358Supplementary analyses
“Rarefaction” analysis
Plot a curve showing how many unique interactions you get with each new video scored:
curve_analysis <- drop_ufos(scored)
curve_analysis <- curve_analysis[complete.cases(curve_analysis),]
curve_analysis <- curve_analysis[,c(4,5,6,10,13)]
# Convert the dates and times into a readable format and sort the dataset by them:
library(chron)
curve_analysis$Date.of.recording <- as.Date(curve_analysis$Date.of.recording, format = "%d %B %y")
# This will makes the times readable by 'chron':
curve_analysis$Time.of.Recording <- paste(curve_analysis$Time.of.Recording, ":00", sep = "")
curve_analysis$Time.of.Recording <- chron(times = curve_analysis$Time.of.Recording)
# Sort by date and time:
curve_analysis <- curve_analysis[order(curve_analysis$Date.of.recording, curve_analysis$Time.of.Recording, decreasing = TRUE),]
species_pairs <- list()
rarefaction <- data.frame(matrix(NA, nrow = nrow(curve_analysis), ncol = 2))
colnames(rarefaction) <- c("Observation", "Unique Pairs")
for(i in 1:nrow(curve_analysis)) {
new_pair <- c(as.character(curve_analysis$Species.A[i]), as.character(curve_analysis$Species.B[i]))
species_pairs <- c(species_pairs, list(new_pair))
rarefaction[i,1] <- i
rarefaction[i,2] <- length(unique(species_pairs))
}
plot(rarefaction$Observation, rarefaction$`Unique Pairs`, pch = 16, xlab = "Observations", ylab = "Unique pairs of interacting species")
# Now order the observations according to scoring time rather than recording time:
curve_analysis2 <- drop_ufos(scored)
curve_analysis2 <- curve_analysis2[complete.cases(curve_analysis2),]
curve_analysis2 <- curve_analysis2[,c(3,10,13)]
curve_analysis2$Date.of.scoring <- as.Date(curve_analysis2$Date.of.scoring, format = "%d %B %Y")
curve_analysis2 <- curve_analysis2[order(curve_analysis2$Date.of.scoring, decreasing = TRUE),]
species_pairs2 <- list()
rarefaction2 <- data.frame(matrix(NA, nrow = nrow(curve_analysis), ncol = 2))
colnames(rarefaction2) <- c("Observation", "Unique Pairs")
for(i in 1:nrow(curve_analysis2)) {
new_pair <- c(as.character(curve_analysis2$Species.A[i]), as.character(curve_analysis2$Species.B[i]))
species_pairs2 <- c(species_pairs2, list(new_pair))
rarefaction2[i,1] <- i
rarefaction2[i,2] <- length(unique(species_pairs2))
}
plot(rarefaction2$Observation, rarefaction2$`Unique Pairs`, pch = 16, xlab = "Observations", ylab = "Unique pairs of interacting species", main = "Contribution of new observations to dataset expansion")
This is the code used to generate Figure S4:
png("/Users/David/Documents/College/2017–18/Winter 2018/FBQ_git/Figure_S1.png", width = 2880, height = 1800, pointsize = 64)
# setEPS()
# postscript("/Users/David/Documents/College/2017–18/Winter 2018/FBQ_git/Figure_S1.eps")
plot(rarefaction2$Observation, rarefaction2$`Unique Pairs`, pch = 16, xlab = "Observations", ylab = "Unique pairs of interacting species", main = "Contribution of new observations to dataset expansion")
dev.off()Breakdown by diet
diet <- read.csv("/Users/David/Documents/College/2017–18/Winter 2018/FBQ_git/diet.csv", sep = ",", strip.white = TRUE)
# Add dietary information for both Species A and Species B
for(i in 1:nrow(both_group)) {
both_group[i, 12] <- diet$Diet[diet$Species == as.character(both_group$Species.A[i])]
both_group[i, 13] <- diet$Diet[diet$Species == as.character(both_group$Species.B[i])]
}
colnames(both_group)[12:13] <- c("Diet.A", "Diet.B")
# Now subsample the test frame to include only those interactions that take place between
# species with the same diet:
test_benth <- both_group[both_group$Diet.A == "Benthic Invertebrate Consumer" &
both_group$Diet.B == "Benthic Invertebrate Consumer",]
nrow(test_benth)## [1] 0test_omni <- both_group[both_group$Diet.A == "Omnivore" & both_group$Diet.B == "Omnivore",]
nrow(test_omni)## [1] 16test_herb <- both_group[both_group$Diet.A == "Herbivore/Detritivore" & both_group$Diet.B == "Herbivore/Detritivore",]
nrow(test_herb)## [1] 4test_coral <- both_group[both_group$Diet.A == "Corallivore" & both_group$Diet.B == "Corallivore",]
nrow(test_coral)## [1] 0# Test the validity of Lanchester's laws for interaction subsets of non-zero length.
# Omnivores:
sum(test_omni$Linear.True)/nrow(test_omni)## [1] 0.8125sum(test_omni$Square.True)/nrow(test_omni)## [1] 0.5625dbinom(sum(test_omni$Linear.True), size = nrow(test_omni), prob = 1/2) ## [1] 0.008544922dbinom(sum(test_omni$Square.True), size = nrow(test_omni), prob = 1/2)## [1] 0.1745605dif3 <- test_omni[test_omni$Linear.Law.Prediction != test_omni$Square.Law.Prediction,]
nrow(dif3)## [1] 4sum(dif3$Linear.True)/nrow(dif3)## [1] 1dbinom(sum(dif3$Linear.True), size = nrow(dif3), prob = 1/2) ## [1] 0.0625# Herbivores/Detritivores:
sum(test_herb$Linear.True)/nrow(test_herb)## [1] 1sum(test_herb$Square.True)/nrow(test_herb)## [1] 0.25dbinom(sum(test_herb$Linear.True), size = nrow(test_herb), prob = 1/2)## [1] 0.0625dbinom(sum(test_herb$Square.True), size = nrow(test_herb), prob = 1/2)## [1] 0.25dif4 <- test_herb[test_herb$Linear.Law.Prediction != test_herb$Square.Law.Prediction,]
nrow(dif4)## [1] 3sum(dif4$Linear.True)/nrow(dif4)## [1] 1dbinom(sum(dif4$Linear.True), size = nrow(dif4), prob = 1/2) ## [1] 0.125# Re-run everything with mass raised to two thirds:
for(i in 1:nrow(both_group2)) {
both_group2[i, 12] <- diet$Diet[diet$Species == as.character(both_group2$Species.A[i])]
both_group2[i, 13] <- diet$Diet[diet$Species == as.character(both_group2$Species.B[i])]
}
colnames(both_group2)[12:13] <- c("Diet.A", "Diet.B")
test_benth2 <- both_group2[both_group2$Diet.A == "Benthic Invertebrate Consumer" &
both_group2$Diet.B == "Benthic Invertebrate Consumer",]
test_omni2 <- both_group2[both_group2$Diet.A == "Omnivore" & both_group2$Diet.B == "Omnivore",]
test_herb2 <- both_group2[both_group2$Diet.A == "Herbivore/Detritivore" & both_group2$Diet.B == "Herbivore/Detritivore",]
sum(test_omni2$Linear.True)/nrow(test_omni2)## [1] 0.5625sum(test_omni2$Square.True)/nrow(test_omni2)## [1] 0.25dbinom(sum(test_omni2$Linear.True), size = nrow(test_omni2), prob = 1/2) ## [1] 0.1745605dbinom(sum(test_omni2$Square.True), size = nrow(test_omni2), prob = 1/2) ## [1] 0.027771dif5 <- test_omni2[test_omni2$Linear.Law.Prediction != test_omni2$Square.Law.Prediction,]
nrow(dif5)## [1] 5sum(dif5$Linear.True)/nrow(dif5)## [1] 1dbinom(sum(dif5$Linear.True), size = nrow(dif5), prob = 1/2) ## [1] 0.03125sum(test_herb2$Linear.True)/nrow(test_herb2)## [1] 0.75sum(test_herb2$Square.True)/nrow(test_herb2)## [1] 0dbinom(sum(test_herb2$Linear.True), size = nrow(test_herb2), prob = 1/2)## [1] 0.25dbinom(sum(test_herb2$Square.True), size = nrow(test_herb2), prob = 1/2)## [1] 0.0625dif6 <- test_herb2[test_herb2$Linear.Law.Prediction != test_herb2$Square.Law.Prediction,]
nrow(dif6)## [1] 3sum(dif6$Linear.True)/nrow(dif6)## [1] 1dbinom(sum(dif6$Linear.True), size = nrow(dif6), prob = 1/2) ## [1] 0.125Analyses not included in the manuscript
These include exploratory analyses as well as applications of previously proposed methods to the fish dataset (e.g., Shelley et al.’s (2004) logistic regression).
Most often interacting species
Subsample the data sheets down to the 10 most often interacting species. This will increase completeness (substantially) and linearity (slightly), but give you fewer data points to work with.
mostoften <- function(sheet, df, n) {
occurrences <- tail(sort(table(unlist(sheet[,c("Species.A", "Species.B")]))), n)
max_n <- as.character(as.data.frame(occurrences)[,1])
max_n_matrix <- subset(df, select = max_n)
max_n_matrix <- subset(max_n_matrix, rownames(max_n_matrix) %in% max_n)
max_n_matrix <- max_n_matrix[order(row.names(max_n_matrix)), order(colnames(max_n_matrix))]
return(max_n_matrix)
}
dommatrix_max10 <- mostoften(sheet4, dommatrix, 10)
scoredmatrix_max10 <- mostoften(scored4, scoredmatrix, 10)
bothmatrix_max10 <- mostoften(both4, bothmatrix, 10)
domooo_max10 <- mostoften(sheet_ooo4, domooo, 10)
scoredooo_max10 <- mostoften(scored_ooo4, scoredooo, 10)
bothooo_max10 <- mostoften(both_ooo4, bothooo, 10)
dommax10 <- convert_to_df(dommatrix_max10)
scoredmax10 <- convert_to_df(scoredmatrix_max10)
bothmax10 <- convert_to_df(bothmatrix_max10)
liveooomax10 <- convert_to_df(domooo_max10)
scoredooomax10 <- convert_to_df(scoredooo_max10)
bothooomax10 <- convert_to_df(bothooo_max10)
paste("The completeness of the live-only 10-species matrix is ", compl(dommax10), "%", sep = "")
paste("The completeness of the video-only 10-species matrix is ", compl(scoredmax10), "%", sep = "")
paste("The completeness of the combined 10-species matrix is ", compl(bothmax10), "%", sep = "")
paste("The completeness of the live-only one-on-one 10-species matrix is ", compl(liveooomax10), "%", sep = "")
paste("The completeness of the video-only one-on-one 10-species matrix is ", compl(scoredooomax10), "%", sep = "")
paste("The completeness of the combined one-on-one 10-species matrix is ", compl(bothooomax10), "%", sep = "")
cbi_max10 <- cbi_hierarchy(dommax10)
cbi_scoredmax10 <- cbi_hierarchy(scoredmax10)
cbi_bothmax10 <- cbi_hierarchy(bothmax10)
cbi_liveooomax10 <- cbi_hierarchy(liveooomax10)
cbi_scoredooomax10 <- cbi_hierarchy(scoredooomax10)
cbi_bothooomax10 <- cbi_hierarchy(bothooomax10)
devries(dommax10, Nperms = 10000, history = FALSE, plot = TRUE)
devries(scoredmax10, Nperms = 10000, history = FALSE, plot = TRUE)
devries(bothmax10, Nperms = 10000, history = FALSE, plot = TRUE)
devries(liveooomax10, Nperms = 10000, history = FALSE, plot = TRUE)
devries(scoredooomax10, Nperms = 10000, history = FALSE, plot = TRUE)
devries(bothooomax10, Nperms = 10000, history = FALSE, plot = TRUE)I&SI rank
A third measure of dominance that is particularly applicable to weakly but significantly linear dominance hierarchies, implemented as described in Schmid & de Vries (2013).
library(compete)
i_and_si_liveooo <- isi13(liveoooall, nTries = 10)
rev(i_and_si_liveooo$best_order) # The order is given from least to most dominant by default
i_and_si_videoooo <- isi13(scoredoooall, nTries = 10)
rev(i_and_si_videoooo$best_order)
i_and_si_bothooo <- isi13(bothoooall, nTries = 10)
rev(i_and_si_bothooo$best_order)
i_and_si <- isi13(dommax10, nTries = 10)
rev(i_and_si$best_order)
i_and_si_scored <- isi13(scoredmax10, nTries = 10)
rev(i_and_si_scored$best_order)
i_and_si_both <- isi13(bothmax10, nTries = 10)
rev(i_and_si_both$best_order)
# Create a data frame with alphabetically ordered species and their corresponding I&SI ranks:
isi_both <- i_and_si_bothooo$best_order
isi_both <- setdiff(isi_both, "Cornetfish")
test_isi <- data.frame(matrix(NA, nrow = length(isi_both), ncol = 2))
colnames(test_isi) <- c("Species", "I_SI_Rank")
test_isi$Species <- sort(isi_both)
test_isi$I_SI_Rank <- match(sort(isi_both), isi_both)Nonlinear models
Do simple nonlinear models (exponential, logarithmic, quadratic) explain the relationship between individual fighting ability (as determined by CBI and the frequency of winning) and body mass?
# Add mass squared to the test frame:
both_test[,9] <- both_test$Mass^2
colnames(both_test)[9] <- "Mass_Squared"
stats <- function(df) {
results <- data.frame(matrix(NA, nrow = 2, ncol = 14))
# Columns represent different statistics associated with each model. Spearman's rho
# and Pearson's coefficient are included for comparison; if the former is significant
# but the latter is not, this shows that the relationship is monotonic but nonlinear.
colnames(results) <- c("rho", "Spear_p-value", "cor", "Pear_p-value", "linAIC", "expRsq", "exp_p-value", "expAIC", "logRsq", "log_p-value", "logAIC", "quadrRsq", "quadr_p-value", "quadrAIC")
# Rows represent dependent variables:
row.names(results) <- c("CBI_Mass", "WinFreq_Mass")
spear_mass_CBI <- cor.test(df$CBI, df$Mass, method = "spearman", alternative = "two.sided")
spear_mass_WinFreq <- cor.test(df$WinFreq, df$Mass, method = "spearman", alternative = "two.sided")
pear_mass_CBI <- cor.test(df$CBI, df$Mass, method = "pearson", alternative = "two.sided")
pear_mass_WinFreq <- cor.test(df$WinFreq, df$Mass, method = "pearson", alternative = "two.sided")
linear_mass_CBI <- lm(df$CBI ~ df$Mass)
expo_mass_CBI <- lm(log(df$CBI) ~ df$Mass)
loga_mass_CBI <- lm(df$CBI ~ log(df$Mass))
quadra_mass_CBI <- lm(df$CBI ~ df$Mass + df$Mass_Squared)
linear_mass_WinFreq <- lm(df$WinFreq ~ df$Mass)
# No exponential model for winning frequency: the minimum winning frequency observed
# in the dataset is 0, and simple log transformation therefore cannot be applied to it.
loga_mass_WinFreq <- lm(df$WinFreq ~ log(df$Mass))
quadra_mass_WinFreq <- lm(df$WinFreq ~ df$Mass + df$Mass_Squared)
results[1,1] <- spear_mass_CBI$estimate; results[1,2] <- spear_mass_CBI$p.value
results[1,3] <- pear_mass_CBI$estimate; results[1,4] <- pear_mass_CBI$p.value
results[1,5] <- AIC(linear_mass_CBI); results[1,6] <- summary(expo_mass_CBI)$r.squared
results[1,7] <- summary(expo_mass_CBI)$coefficients[,4][[2]]
results[1,8] <- AIC(expo_mass_CBI); results[1,9] <- summary(loga_mass_CBI)$r.squared
results[1,10] <- summary(loga_mass_CBI)$coefficients[,4][[2]]
results[1,11] <- AIC(loga_mass_CBI); results[1,12] <- summary(quadra_mass_CBI)$r.squared
results[1,13] <- summary(quadra_mass_CBI)$coefficients[,4][[2]]
results[1,14] <- AIC(quadra_mass_CBI)
results[2,1] <- spear_mass_WinFreq$estimate; results[2,2] <- spear_mass_WinFreq$p.value
results[2,3] <- pear_mass_WinFreq$estimate; results[2,4] <- pear_mass_WinFreq$p.value
results[2,5] <- AIC(linear_mass_WinFreq)
results[2,9] <- summary(loga_mass_WinFreq)$r.squared
results[2,10] <- summary(loga_mass_WinFreq)$coefficients[,4][[2]]
results[2,11] <- AIC(loga_mass_WinFreq)
results[2,12] <- summary(quadra_mass_WinFreq)$r.squared
results[2,13] <- summary(quadra_mass_WinFreq)$coefficients[,4][[2]]
results[2,14] <- AIC(quadra_mass_WinFreq)
return(results)
}
both_test_results <- stats(both_test)
# Adjust the p-values for multiple comparisons:
library(stats)
pvalues <- unlist(both_test_results[,c(2,4,7,10,13)])
p.adjust(pvalues, "bonferroni", sum(!is.na(pvalues)))Logistic regression
The code below is an implementation of the logistic regression analysis described by Shelley et al. (2004). Given the paucity of data, no pairwise comparisons were made; each of the selected species was compared against all other species it interacted with. Unlike Shelley et al. (2004), we regress against both group size and group size squared.
# Extract those species that participated in an interaction with more than one individual
# at least a given number of times:
atleastoncemorethan <- function(df, n) {
output <- vector()
species <- union(unique(df$Species.A), unique(df$Species.B))
for (i in species) {
a.frame <- df[df$Species.A == i,]
b.frame <- df[df$Species.B == i,]
if (sum(!is.na(a.frame$Group.A) & a.frame$Group.A > 1) + sum(!is.na(b.frame$Group.B) & b.frame$Group.B > 1) >= n) {
output <- c(output, as.character(i))
}
}
return(unique(output))
}
# Some ad hoc data clean-up before applying the atleastoncemorethan() function:
scored4$Group.B <- as.numeric(as.character(scored4$Group.B))
both4$Group.B <- as.numeric(both4$Group.B)
# Include only those species that had a group size greater than 1 at least twice. This
# ensures that there is variability in the predictor variable.
live_logres <- sheet4[sheet4$Species.A %in% atleastoncemorethan(sheet4, 2) &
sheet4$Species.B %in% atleastoncemorethan(sheet4, 2),]
video_logres <- scored4[scored4$Species.A %in% atleastoncemorethan(scored4, 2) &
scored4$Species.B %in% atleastoncemorethan(scored4, 2),]
both_logres <- both4[both4$Species.A %in% atleastoncemorethan(both4, 2) &
both4$Species.B %in% atleastoncemorethan(both4, 2),]
# Convert from wide to long format: every two successive rows correspond to one
# interaction, with the winner coded as "1" and the loser coded as "0".
wide_to_long <- function(df) {
longframe <- data.frame(matrix(NA, ncol = 3, nrow = 2*nrow(df)))
colnames(longframe) <- c("Species", "Group_Size", "Is_Winner")
for (i in 1:nrow(df)) {
longframe$Species[2*i-1] <- as.character(df$Species.A[i])
longframe$Species[2*i] <- as.character(df$Species.B[i])
longframe$Group_Size[2*i-1] <- df$Group.A[i]
longframe$Group_Size[2*i] <- df$Group.B[i]
if (identical(as.character(df$Species.A[i]), as.character(df$Winner[i])) == TRUE) {
longframe$Is_Winner[2*i-1] <- 1
longframe$Is_Winner[2*i] <- 0
}
if (identical(as.character(df$Species.B[i]), as.character(df$Winner[i])) == TRUE) {
longframe$Is_Winner[2*i-1] <- 0
longframe$Is_Winner[2*i] <- 1
}
}
return(longframe)
}
# Make sure winner is always either A or B, and that you have group sizes for all rows:
live_long <- wide_to_long(live_logres)
nrow(live_long) == nrow(live_long[!is.na(live_long$Is_Winner),])
nrow(live_long) == nrow(live_long[!is.na(live_long$Group_Size),])
video_long <- wide_to_long(video_logres)
nrow(video_long) == nrow(video_long[!is.na(video_long$Is_Winner),])
nrow(video_long) == nrow(video_long[!is.na(video_long$Group_Size),])
both_long <- wide_to_long(both_logres)
nrow(both_long) == nrow(both_long[!is.na(both_long$Is_Winner),])
nrow(both_long) == nrow(both_long[!is.na(both_long$Group_Size),])
# Add group size squared
square_it <- function(df) {
a <- ncol(df)
df[,a+1] <- df$Group_Size^2
colnames(df)[a+1] <- "Group_Size_Squared"
return(df)
}
live_long <- square_it(live_long)
video_long <- square_it(video_long)
both_long <- square_it(both_long)
# Extract those species that won and lost at least 1 interaction. This ensures that there
# is variability in the dependent variable:
won_once <- function(df) {
output <- vector()
species <- unique(df$Species)
for(i in species) {
subset.frame <- df[df$Species == i,]
if(0 %in% subset.frame$Is_Winner & 1 %in% subset.frame$Is_Winner) {
output <- c(output, i)
}
}
return(df[df$Species %in% output,])
}
# This function creates a data frame with logistic regression results. The regression is
# run against group size if x is FALSE and against group size squared if x is TRUE.
library(pscl)
species_logistic <- function(df, x = TRUE) {
j <- 1
species <- unique(df$Species)
results <- data.frame(matrix(NA, nrow = length(species), ncol = 6))
colnames(results) <- c("Species", "Coefficient", "p-value", "McFaddenpseudoRsq", "AIC", "n")
for(i in species) {
if (x == FALSE) {
model <- glm(df$Is_Winner[df$Species == i] ~ df$Group_Size[df$Species == i])
} else {
model <- glm(df$Is_Winner[df$Species == i] ~ df$Group_Size_Squared[df$Species == i] + df$Group_Size[df$Species == i])
}
results[j,1] <- i
results[j,2] <- model$coefficients[2]
results[j,3] <- coef(summary(model))[,4][2]
results[j,4] <- pR2(model)[4]
results[j,5] <- AIC(model)
results[j,6] <- nrow(df[df$Species == i,])
j <- j + 1
}
return(results)
}
# Linear group size
species_logistic(won_once(live_long), FALSE)
species_logistic(won_once(video_long), FALSE)
species_logistic(won_once(both_long), FALSE)
# Group size squared
species_logistic(won_once(live_long), TRUE)
species_logistic(won_once(video_long), TRUE)
species_logistic(won_once(both_long), TRUE)Chi-squared
This is an alternative way of testing Lanchester’s laws. Instead of applying the binomial test to a binary variable (prediction: correct, incorrect), we apply \(\chi^2\) to a 4-entry contingency table (correctly predicted win, incorrectly predicted win, correctly predicted loss, incorrectly predicted loss):
populate_cont_table <- function(df, law) {
# Create a contingency table whose rows represent (1) Wins, (2) Losses predicted
# by the linear law and whose columns represent observed (1) Wins and (2) Losses
cont_table <- matrix(0, nrow = 2, ncol = 2)
# Focus solely on Species A. This is arbitrary; we could instead focus on B by
# interchanging individual rows and individual columns.
for(i in 1:nrow(df)) {
# A won:
if (identical(as.character(df$Species.A[i]), as.character(df$Winner[i])) == TRUE) {
# A was predicted to win:
if (identical(as.character(df$Species.A[i]), as.character(df[i, paste(law, ".Law.Prediction", sep = "")])) == TRUE) {
inc(cont_table[1,1])
# A was predicted to lose:
} else {
inc(cont_table[2,1])
}
# A lost:
} else {
# A was predicted to win:
if (identical(as.character(df$Species.A[i]), as.character(df[i, paste(law, ".Law.Prediction", sep = "")])) == TRUE) {
inc(cont_table[1,2])
# A was predicted to lose:
} else {
inc(cont_table[2,2])
}
}
}
return(cont_table)
}
chisq.test(populate_cont_table(both_group, "Linear"))
chisq.test(populate_cont_table(both_group, "Square"))
chisq.test(populate_cont_table(both_group2, "Linear"))
chisq.test(populate_cont_table(both_group2, "Square"))
# To distinguish between the linear and square laws, we run another chi-squared test.
# This time, the contingency table has the following structure: [1,1] correctly
# predicted by both models, [1,2] correctly predicted by the linear law but not the
# square law, [2,1] correctly predicted by the square law but not the linear law,
# [2,2] predicted incorrectly by both laws.
law_vs_law <- function(df) {
cont_table <- matrix(0, nrow = 2, ncol = 2)
for(i in 1:nrow(df)) {
if (df$Linear.True[i] == 1 & df$Square.True[i] == 1) {
inc(cont_table[1,1])
} else if (df$Linear.True[i] == 1 & df$Square.True[i] == 0) {
inc(cont_table[1,2])
} else if (df$Linear.True[i] == 0 & df$Square.True[i] == 1) {
inc(cont_table[2,1])
} else {
inc(cont_table[2,2])
}
}
return(cont_table)
}
laws_mass <- law_vs_law(both_group)
laws_mass2_3 <- law_vs_law(both_group2)
chisq.test(laws_mass)
chisq.test(laws_mass2_3)Response to Behav. Ecol. reviews
Fish size correlation
We discovered a discrepancy between the identification of the longnose butterflyfish in our fish size data (FishSize.csv) and Table 1 of our manuscript: the same species was identified as Forcipiger longirostris in the former and as the closely related F. flavissimus in the latter. Based on the photograph we took of the species, we concluded that the latter identification was correct, as the lower half of its eye was silver rather than black (http://reefs.com/2015/12/31/long-nosed-butterflies-part-2-forcipiger/). We therefore replaced the maximum attained size reported for F. longirostris by Nelson (2005) with that reported for F. flavissimus. Below, we recalculate the correlation between the reported attained sizes and the observed sizes:
# Maximum attained sizes according to Randall (2005):
sizes <- read.csv("/Users/David/Documents/College/2017–18/Winter 2018/FBQ_git/FishSize.csv", sep = ";", header = TRUE)
colnames(sizes)[1] <- "Species"
booksizes <- sizes$Attain[c(1:3,5:8,10:24)]
# Distributions of sizes observed in the field:
fishsize <- read.csv("/Users/David/Documents/College/2017–18/Winter 2018/FBQ_git/FieldFishSize.csv", sep = ",", strip.white = TRUE)
# Calculate the mean, median, and maximum of the observed size distributions:
meansizes <- colMeans(fishsize[,c(2:4,6:9,11:25)], na.rm = TRUE)
colMedians <- function(matrix) {
return(apply(matrix, 2, median, na.rm = TRUE))
}
colMax <- function(matrix) {
return(apply(matrix, 2, max, na.rm = TRUE))
}
mediansizes <- colMedians(fishsize[,c(2:4,6:9,11:25)])
maxsizes <- colMax(fishsize[,c(2:4,6:9,11:25)])
# Do book sizes correlate with live observation sizes?
g <- cor.test(booksizes, mediansizes, method = "spearman", alternative = "two.sided")## Warning in cor.test.default(booksizes, mediansizes, method = "spearman", :
## Cannot compute exact p-value with tiespaste("correlation: ", signif(g$estimate[[1]],5), ", p-value: ", g$p.value, sep = "")## [1] "correlation: 0.71489, p-value: 0.000184888348145873"h <- cor.test(booksizes, meansizes, method = "spearman", alternative = "two.sided")## Warning in cor.test.default(booksizes, meansizes, method = "spearman",
## alternative = "two.sided"): Cannot compute exact p-value with tiespaste("correlation: ", signif(h$estimate[[1]],5), ", p-value: ", h$p.value, sep = "")## [1] "correlation: 0.76279, p-value: 3.66099985642471e-05"i <- cor.test(booksizes, maxsizes, method = "spearman", alternative = "two.sided")## Warning in cor.test.default(booksizes, maxsizes, method = "spearman",
## alternative = "two.sided"): Cannot compute exact p-value with tiespaste("correlation: ", signif(i$estimate[[1]],5), ", p-value: ", i$p.value, sep = "")## [1] "correlation: 0.81934, p-value: 3.09742996782284e-06"j <- cor.test(booksizes, mediansizes, method = "pearson", alternative = "two.sided")
paste("correlation: ", signif(j$estimate[[1]],5), ", p-value: ", j$p.value, sep = "")## [1] "correlation: 0.83931, p-value: 1.04831251950256e-06"k <- cor.test(booksizes, meansizes, method = "pearson", alternative = "two.sided")
paste("correlation: ", signif(k$estimate[[1]],5), ", p-value: ", k$p.value, sep = "")## [1] "correlation: 0.85355, p-value: 4.41082818979503e-07"l <- cor.test(booksizes, maxsizes, method = "pearson", alternative = "two.sided")
paste("correlation: ", signif(l$estimate[[1]],5), ", p-value: ", l$p.value, sep = "")## [1] "correlation: 0.83313, p-value: 1.48832376900543e-06"CBI Calculation
Going back to Clutton-Brock et al. (1979), we discovered that our original implementation of the CBI-calculating function had an error. When counting the species beaten by those species that were dominated by the subject, we need to exclude the subject itself; the same is true when tallying the species that dominated those species which dominated the subject. This was not done in the original code, but it is done below:
# Remember:
# The A column represents the species that A dominates (A wins)
# The A row represents the species that dominate A (A losses)
# (B + b + 1)/(L + l + 1)
cbi_new <- function(df, species) {
# number of entries in column 'species' that are neither 0 nor NA:
B <- sum(df[,species] != 0, na.rm = TRUE)
# species corresponding to those entries:
dominates <- rownames(df)[which(df[,species] != 0)]
# number of entries in the columns corresponding to those species that are neither 0 nor NA:
if (length(dominates) > 0) {
# exclude 'species' itself from consideration
b <- sum(df[rownames(df) != species, colnames(df) %in% species] != 0, na.rm = TRUE)
} else {
b <- 0
}
# number of entries in row B that are neither 0 nor NA:
L <- sum(df[species,], na.rm = TRUE)
# species corresponding to those entries:
dominated <- colnames(df)[which(df[species,] != 0)]
# number of entries in the rows corresponding to those species that are neither 0 nor NA:
if (length(dominated) > 0) {
# exclude 'species' itself from consideration
l <- sum(df[rownames(df) %in% dominated, colnames(df) != species] != 0, na.rm = TRUE)
} else {
l <- 0
}
return((B + b + 1)/(L + l + 1))
}
# Apply it to the full matrix:
cbi_hierarchy_new <- function(df) {
cbi_matrix <- data.frame(matrix(NA, nrow = nrow(df), ncol = 2))
colnames(cbi_matrix) <- c("Species", "CBI")
for (i in 1:nrow(df)) {
cbi_matrix$Species[i] <- rownames(df)[i]
cbi_matrix$CBI[i] <- cbi_new(df, rownames(df)[i])
}
cbi_matrix <- cbi_matrix[order(cbi_matrix$CBI, decreasing = TRUE),]
return(cbi_matrix)
}Phylogenetic GLS
First, we need to drop the “Juvenile Parrotfish” category from the analysis, since it is a mixture of multiple species:
nojp_domall <- exclude(nocornet_domall, "Juvenile Parrotfish")
nojp_scoredall <- exclude(scoredall, "Juvenile Parrotfish")
nojp_bothall <- exclude(nocornet_bothall, "Juvenile Parrotfish")
nojp_live <- exclude(nocornet_live, "Juvenile Parrotfish")
nojp_scored <- exclude(scoredoooall, "Juvenile Parrotfish")
nojp_both <- exclude(nocornet_both, "Juvenile Parrotfish")
# Recalculate CBI:
cbi_nojp_domall <- cbi_hierarchy_new(nojp_domall)
cbi_nojp_scoredall <- cbi_hierarchy_new(nojp_scoredall)
cbi_nojp_bothall <- cbi_hierarchy_new(nojp_bothall)
cbi_nojp_ooolive <- cbi_hierarchy_new(nojp_live)
cbi_nojp_ooovideo <- cbi_hierarchy_new(nojp_scored)
cbi_nojp_oooboth <- cbi_hierarchy_new(nojp_both)
# Recalculate the frequency of winning and combine the results with the CBI hierarchy
# and mass-length data into a single data frame for further testing. We will again focus
# on one-on-one data only from now on.
# This is a generalized version of the no_surgeon() function defined above:
drop_taxa <- function(df, taxa_to_drop) {
return(df[!df$Species.A %in% taxa_to_drop & !df$Species.B %in% taxa_to_drop,])
}
sheet_ooo5 <- drop_taxa(sheet_ooo4, c("Cornetfish", "Juvenile Parrotfish"))
scored_ooo5 <- drop_taxa(scored_ooo4, c("Cornetfish", "Juvenile Parrotfish"))
both_ooo5 <- drop_taxa(both_ooo4, c("Cornetfish", "Juvenile Parrotfish"))
sheet_phylo_test <- merge(cbi_nojp_ooolive, win_freq(sheet_ooo5), by = "Species")
sheet_phylo_test <- merge(sheet_phylo_test, lengthmass, by = "Species")
# We need to provide row names to ensure that each row of the data frame is linked to the
# correct tree tip:
rownames(sheet_phylo_test) <- sheet_phylo_test$Species
scored_phylo_test <- merge(cbi_nojp_ooovideo, win_freq(scored_ooo5), by = "Species")
scored_phylo_test <- merge(scored_phylo_test, lengthmass, by = "Species")
rownames(scored_phylo_test) <- scored_phylo_test$Species
both_phylo_test <- merge(cbi_nojp_oooboth, win_freq(both_ooo5), by = "Species")
both_phylo_test <- merge(both_phylo_test, lengthmass, by = "Species")
rownames(both_phylo_test) <- both_phylo_test$Species
# Finally, repeat all of the above without the spotted porcupinefish (we will want to see
# how its exclusion as an outlier is going to change the results). This species was never
# observed on video, so there is no need to remove it from the video-only data. Note that
# the exclude() function is not safeguarded against attempts to remove taxa that are not
# present in the data frame to begin with.
nojpsp_domall <- exclude(nojp_domall, "Spotted Porcupinefish")
nojpsp_bothall <- exclude(nojp_bothall, "Spotted Porcupinefish")
nojpsp_live <- exclude(nojp_live, "Spotted Porcupinefish")
nojpsp_both <- exclude(nojp_both, "Spotted Porcupinefish")
cbi_nojpsp_domall <- cbi_hierarchy_new(nojpsp_domall)
cbi_nojpsp_bothall <- cbi_hierarchy_new(nojpsp_bothall)
cbi_nojpsp_ooolive <- cbi_hierarchy_new(nojpsp_live)
cbi_nojpsp_oooboth <- cbi_hierarchy_new(nojpsp_both)
sheet_ooo6 <- drop_taxa(sheet_ooo4, c("Cornetfish", "Juvenile Parrotfish", "Spotted Porcupinefish"))
both_ooo6 <- drop_taxa(both_ooo4, c("Cornetfish", "Juvenile Parrotfish", "Spotted Porcupinefish"))
sheet_phylo_test_nosp <- merge(cbi_nojpsp_ooolive, win_freq(sheet_ooo6), by = "Species")
sheet_phylo_test_nosp <- merge(sheet_phylo_test_nosp, lengthmass, by = "Species")
rownames(sheet_phylo_test_nosp) <- sheet_phylo_test_nosp$Species
both_phylo_test_nosp <- merge(cbi_nojpsp_oooboth, win_freq(both_ooo6), by = "Species")
both_phylo_test_nosp <- merge(both_phylo_test_nosp, lengthmass, by = "Species")
rownames(both_phylo_test_nosp) <- both_phylo_test_nosp$SpeciesRecalculate the properties of the one-on-one matrices created by the exclusion of juvenile parrotfish:
paste("The completeness of the live-only one-on-one matrix after the exclusion of juvenile parrotfish is ", compl(nojp_live), "%", sep = "")## [1] "The completeness of the live-only one-on-one matrix after the exclusion of juvenile parrotfish is 35.7142857142857%"paste("The completeness of the video-only one-on-one matrix after the exclusion of juvenile parrotfish is ", compl(nojp_scored), "%", sep = "")## [1] "The completeness of the video-only one-on-one matrix after the exclusion of juvenile parrotfish is 26.6666666666667%"paste("The completeness of the combined one-on-one matrix after the exclusion of juvenile parrotfish is ", compl(nojp_both), "%", sep = "")## [1] "The completeness of the combined one-on-one matrix after the exclusion of juvenile parrotfish is 34.6320346320346%"devries(nojp_live, Nperms = 10000, history = FALSE, plot = TRUE)## h-modified = 0.3125917
## p-value= 0.0042
devries(nojp_scored, Nperms = 10000, history = FALSE, plot = TRUE)## h-modified = 0.3038453
## p-value= 0.07
devries(nojp_both, Nperms = 10000, history = FALSE, plot = TRUE)## h-modified = 0.2824597
## p-value= 0.0062
Now we can run phylogenetic regression:
library(nlme)
library(phytools)## Loading required package: ape## Loading required package: mapsfishtree <- read.tree("/Users/David/Documents/College/2017–18/Winter 2018/FBQ paper/Behav Ecol revisions/small_tree.tre")
# We will have to match the scientific names in the tip labels against the common names
# we use for the remaining data:
species_table <- read.csv("/Users/David/Documents/College/2017–18/Winter 2018/FBQ_git/FishSize.csv", sep = ";", stringsAsFactors = FALSE)
fishtree$tip.label <- species_table$Common.Name[match(unlist(lapply(strsplit(fishtree$tip.label, "_"), function(x) paste(x[1], x[2], sep = " "))), species_table$Scientific.Name)]
# Make a new tree that does not include the spotted porcupinefish:
fishtree_nosp <- drop.tip(fishtree, "Spotted Porcupinefish")
# Run phylogenetic GLS on all parameter combinations (mass versus mass raised to two thirds, CBI vs. frequency of winning, spotted porcupinefish vs. no spotted porcupinefish):
pgls_CBI_Mass <- gls(CBI ~ Mass, correlation = corPagel(1, phy = fishtree), data = both_phylo_test, method = "ML")
summary(pgls_CBI_Mass)## Generalized least squares fit by maximum likelihood
## Model: CBI ~ Mass
## Data: both_phylo_test
## AIC BIC logLik
## 42.96499 47.32916 -17.4825
##
## Correlation Structure: corPagel
## Formula: ~1
## Parameter estimate(s):
## lambda
## 1.026069
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) 0.6181935 0.3625948 1.7049155 0.1037
## Mass -0.0002343 0.0006324 -0.3705309 0.7149
##
## Correlation:
## (Intr)
## Mass -0.341
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -0.7446376 -0.6191146 -0.3996205 -0.2396558 3.0142764
##
## Residual standard error: 0.8129451
## Degrees of freedom: 22 total; 20 residualpgls_CBI_Mass2_3 <- gls(CBI ~ Mass2_3, correlation = corPagel(1, phy = fishtree), data = both_phylo_test, method = "ML")
summary(pgls_CBI_Mass2_3)## Generalized least squares fit by maximum likelihood
## Model: CBI ~ Mass2_3
## Data: both_phylo_test
## AIC BIC logLik
## 43.09373 47.4579 -17.54687
##
## Correlation Structure: corPagel
## Formula: ~1
## Parameter estimate(s):
## lambda
## 1.024889
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) 0.5973099 0.3862164 1.5465678 0.1376
## Mass2_3 -0.0008940 0.0065095 -0.1373374 0.8921
##
## Correlation:
## (Intr)
## Mass2_3 -0.47
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -0.7163492 -0.5910616 -0.4011444 -0.2178379 3.0026720
##
## Residual standard error: 0.8133159
## Degrees of freedom: 22 total; 20 residualpgls_WinFreq_Mass <- gls(WinFreq ~ Mass, correlation = corPagel(1, phy = fishtree), data = both_phylo_test, method = "ML")
summary(pgls_WinFreq_Mass)## Generalized least squares fit by maximum likelihood
## Model: WinFreq ~ Mass
## Data: both_phylo_test
## AIC BIC logLik
## 9.610968 13.97514 -0.8054841
##
## Correlation Structure: corPagel
## Formula: ~1
## Parameter estimate(s):
## lambda
## 0.7608329
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) 0.3867176 0.12406998 3.1169316 0.0054
## Mass 0.0001547 0.00024972 0.6193854 0.5427
##
## Correlation:
## (Intr)
## Mass -0.389
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -1.3855278 -0.9836204 -0.3401884 0.1149432 1.9053666
##
## Residual standard error: 0.2980929
## Degrees of freedom: 22 total; 20 residualpgls_WinFreq_Mass2_3 <- gls(WinFreq ~ Mass2_3, correlation = corPagel(1, phy = fishtree), data = both_phylo_test, method = "ML")
summary(pgls_WinFreq_Mass2_3)## Generalized least squares fit by maximum likelihood
## Model: WinFreq ~ Mass2_3
## Data: both_phylo_test
## AIC BIC logLik
## 9.341281 13.70545 -0.6706403
##
## Correlation Structure: corPagel
## Formula: ~1
## Parameter estimate(s):
## lambda
## 0.7418624
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) 0.3558832 0.13389644 2.6578989 0.0151
## Mass2_3 0.0021802 0.00266542 0.8179602 0.4230
##
## Correlation:
## (Intr)
## Mass2_3 -0.551
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -1.4406599 -0.9320170 -0.4012155 0.1434963 1.9720109
##
## Residual standard error: 0.2934742
## Degrees of freedom: 22 total; 20 residualpgls_CBI_Mass_nosp <- gls(CBI ~ Mass, correlation = corPagel(1, phy = fishtree_nosp), data = both_phylo_test_nosp, method = "ML")
summary(pgls_CBI_Mass_nosp)## Generalized least squares fit by maximum likelihood
## Model: CBI ~ Mass
## Data: both_phylo_test_nosp
## AIC BIC logLik
## 38.77226 42.95035 -15.38613
##
## Correlation Structure: corPagel
## Formula: ~1
## Parameter estimate(s):
## lambda
## 0.9780776
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) 0.4131129 0.3416725 1.209090 0.2415
## Mass 0.0017121 0.0013866 1.234727 0.2320
##
## Correlation:
## (Intr)
## Mass -0.476
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -0.9737380 -0.6492455 -0.4673056 -0.2472223 2.9127563
##
## Residual standard error: 0.7159469
## Degrees of freedom: 21 total; 19 residualpgls_CBI_Mass2_3_nosp <- gls(CBI ~ Mass2_3, correlation = corPagel(1, phy = fishtree_nosp), data = both_phylo_test_nosp, method = "ML")
summary(pgls_CBI_Mass2_3_nosp)## Generalized least squares fit by maximum likelihood
## Model: CBI ~ Mass2_3
## Data: both_phylo_test_nosp
## AIC BIC logLik
## 39.20199 43.38008 -15.60099
##
## Correlation Structure: corPagel
## Formula: ~1
## Parameter estimate(s):
## lambda
## 0.9938416
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) 0.4139249 0.3738005 1.1073418 0.2820
## Mass2_3 0.0093035 0.0095417 0.9750351 0.3418
##
## Correlation:
## (Intr)
## Mass2_3 -0.551
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -0.8003119 -0.6502855 -0.5034040 -0.2894925 2.9422145
##
## Residual standard error: 0.7394883
## Degrees of freedom: 21 total; 19 residualpgls_WinFreq_Mass_nosp <- gls(WinFreq ~ Mass, correlation = corPagel(1, phy = fishtree_nosp), data = both_phylo_test_nosp, method = "ML")
summary(pgls_WinFreq_Mass_nosp)## Generalized least squares fit by maximum likelihood
## Model: WinFreq ~ Mass
## Data: both_phylo_test_nosp
## AIC BIC logLik
## 9.614419 13.79251 -0.8072095
##
## Correlation Structure: corPagel
## Formula: ~1
## Parameter estimate(s):
## lambda
## 0.7513216
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) 0.3365943 0.13968518 2.4096639 0.0263
## Mass 0.0006368 0.00066595 0.9561862 0.3510
##
## Correlation:
## (Intr)
## Mass -0.554
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -1.4891060 -0.9730536 -0.5088269 0.0878157 2.0020512
##
## Residual standard error: 0.2987457
## Degrees of freedom: 21 total; 19 residualpgls_WinFreq_Mass2_3_nosp <- gls(WinFreq ~ Mass2_3, correlation = corPagel(1, phy = fishtree_nosp), data = both_phylo_test_nosp, method = "ML")
summary(pgls_WinFreq_Mass2_3_nosp)## Generalized least squares fit by maximum likelihood
## Model: WinFreq ~ Mass2_3
## Data: both_phylo_test_nosp
## AIC BIC logLik
## 9.600858 13.77895 -0.8004288
##
## Correlation Structure: corPagel
## Formula: ~1
## Parameter estimate(s):
## lambda
## 0.7461108
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) 0.31370935 0.15285260 2.0523652 0.0542
## Mass2_3 0.00450538 0.00465324 0.9682252 0.3451
##
## Correlation:
## (Intr)
## Mass2_3 -0.653
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -1.5174832 -1.0098362 -0.5279994 0.0544028 1.9944100
##
## Residual standard error: 0.2978518
## Degrees of freedom: 21 total; 19 residualWe can see that for all models involving CBI, we get very weak relationships (characterized by nearly horizontal regression lines) that are nevertheless extremely significant. A more careful examination reveals that this is because Pagel’s \(\lambda\) estimated for these models exceeds 1 (moreover, it is identical in all the four cases), even though it is only defined between 0 and 1. The cause of this is unclear (but likely related to the shape of the tree, especially to the presence of terminal branches longer than expected under the Yule process), as is the interpretability of the resulting estimate. Based on Liam Revell’s remarks in the first of the two links given below, we turn to the caper package to fit a phylogenetic GLS under the constraint that \(\lambda \leq 1\):
http://blog.phytools.org/2016/01/phylogenetic-regression-when-estimated.html http://grokbase.com/t/r/r-sig-phylo/113s5mhjgg/pagels-lambda-greater-than-one
if(!require("caper")) {install.packages("caper")}## Loading required package: caper## Loading required package: MASS## Loading required package: mvtnormlibrary(caper)
both_comparative <- comparative.data(data = both_phylo_test, phy = fishtree, names.col = "Species", vcv = TRUE, warn.dropped = TRUE)
both_comparative_nosp <- comparative.data(data = both_phylo_test_nosp, phy = fishtree_nosp, names.col = "Species", vcv = TRUE, warn.dropped = TRUE)
pgls_CBI_Mass <- pgls(formula = CBI ~ Mass, data = both_comparative, lambda = "ML")
summary(pgls_CBI_Mass)##
## Call:
## pgls(formula = CBI ~ Mass, data = both_comparative, lambda = "ML")
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.14667 -0.02814 0.01556 0.05592 0.08199
##
## Branch length transformations:
##
## kappa [Fix] : 1.000
## lambda [ ML] : 1.000
## lower bound : 0.000, p = 0.010507
## upper bound : 1.000, p = 1
## 95.0% CI : (0.524, NA)
## delta [Fix] : 1.000
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.61306817 0.34503191 1.7768 0.09081 .
## Mass -0.00020831 0.00061615 -0.3381 0.73882
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.07567 on 20 degrees of freedom
## Multiple R-squared: 0.005683, Adjusted R-squared: -0.04403
## F-statistic: 0.1143 on 1 and 20 DF, p-value: 0.7388pgls_CBI_Mass2_3 <- pgls(formula = CBI ~ Mass2_3, data = both_comparative, lambda = "ML")
summary(pgls_CBI_Mass2_3)##
## Call:
## pgls(formula = CBI ~ Mass2_3, data = both_comparative, lambda = "ML")
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.14586 -0.02786 0.01345 0.05485 0.09016
##
## Branch length transformations:
##
## kappa [Fix] : 1.000
## lambda [ ML] : 1.000
## lower bound : 0.000, p = 0.016049
## upper bound : 1.000, p = 1
## 95.0% CI : (0.453, NA)
## delta [Fix] : 1.000
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.58731317 0.37068156 1.5844 0.1288
## Mass2_3 -0.00053638 0.00644425 -0.0832 0.9345
##
## Residual standard error: 0.07587 on 20 degrees of freedom
## Multiple R-squared: 0.0003463, Adjusted R-squared: -0.04964
## F-statistic: 0.006928 on 1 and 20 DF, p-value: 0.9345pgls_CBI_Mass_nosp <- pgls(formula = CBI ~ Mass, data = both_comparative_nosp, lambda = "ML")
summary(pgls_CBI_Mass_nosp)##
## Call:
## pgls(formula = CBI ~ Mass, data = both_comparative_nosp, lambda = "ML")
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.14830 -0.04944 -0.01476 0.01939 0.18616
##
## Branch length transformations:
##
## kappa [Fix] : 1.000
## lambda [ ML] : 0.978
## lower bound : 0.000, p = 0.0081072
## upper bound : 1.000, p = 0.83889
## 95.0% CI : (0.458, NA)
## delta [Fix] : 1.000
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.4131104 0.3416712 1.2091 0.2415
## Mass 0.0017121 0.0013866 1.2347 0.2320
##
## Residual standard error: 0.06986 on 19 degrees of freedom
## Multiple R-squared: 0.07428, Adjusted R-squared: 0.02556
## F-statistic: 1.525 on 1 and 19 DF, p-value: 0.232pgls_CBI_Mass2_3_nosp <- pgls(formula = CBI ~ Mass2_3, data = both_comparative_nosp, lambda = "ML")
summary(pgls_CBI_Mass2_3_nosp)##
## Call:
## pgls(formula = CBI ~ Mass2_3, data = both_comparative_nosp, lambda = "ML")
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.151213 -0.018977 0.009497 0.064289 0.139694
##
## Branch length transformations:
##
## kappa [Fix] : 1.000
## lambda [ ML] : 0.994
## lower bound : 0.000, p = 0.0090652
## upper bound : 1.000, p = 0.95565
## 95.0% CI : (0.461, NA)
## delta [Fix] : 1.000
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.4139220 0.3737994 1.1073 0.2820
## Mass2_3 0.0093036 0.0095417 0.9750 0.3418
##
## Residual standard error: 0.07216 on 19 degrees of freedom
## Multiple R-squared: 0.04765, Adjusted R-squared: -0.00247
## F-statistic: 0.9507 on 1 and 19 DF, p-value: 0.3418pgls_WinFreq_Mass <- pgls(formula = WinFreq ~ Mass, data = both_comparative, lambda = "ML")
summary(pgls_WinFreq_Mass)##
## Call:
## pgls(formula = WinFreq ~ Mass, data = both_comparative, lambda = "ML")
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.067829 -0.016253 0.002191 0.017853 0.053835
##
## Branch length transformations:
##
## kappa [Fix] : 1.000
## lambda [ ML] : 0.761
## lower bound : 0.000, p = 0.026067
## upper bound : 1.000, p = 0.060639
## 95.0% CI : (0.083, NA)
## delta [Fix] : 1.000
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.38671755 0.12406983 3.1169 0.005432 **
## Mass 0.00015467 0.00024972 0.6194 0.542654
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.02902 on 20 degrees of freedom
## Multiple R-squared: 0.01882, Adjusted R-squared: -0.03024
## F-statistic: 0.3836 on 1 and 20 DF, p-value: 0.5427pgls_WinFreq_Mass2_3 <- pgls(formula = WinFreq ~ Mass2_3, data = both_comparative, lambda = "ML")
summary(pgls_WinFreq_Mass2_3)##
## Call:
## pgls(formula = WinFreq ~ Mass2_3, data = both_comparative, lambda = "ML")
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.056556 -0.021804 -0.006782 0.014359 0.038175
##
## Branch length transformations:
##
## kappa [Fix] : 1.000
## lambda [ ML] : 0.742
## lower bound : 0.000, p = 0.028143
## upper bound : 1.000, p = 0.052267
## 95.0% CI : (0.071, NA)
## delta [Fix] : 1.000
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.3558831 0.1338963 2.6579 0.0151 *
## Mass2_3 0.0021802 0.0026654 0.8180 0.4230
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.02857 on 20 degrees of freedom
## Multiple R-squared: 0.03237, Adjusted R-squared: -0.01601
## F-statistic: 0.6691 on 1 and 20 DF, p-value: 0.423pgls_WinFreq_Mass_nosp <- pgls(formula = WinFreq ~ Mass, data = both_comparative_nosp, lambda = "ML")
summary(pgls_WinFreq_Mass_nosp)##
## Call:
## pgls(formula = WinFreq ~ Mass, data = both_comparative_nosp,
## lambda = "ML")
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.063430 -0.018323 -0.003256 0.012007 0.054953
##
## Branch length transformations:
##
## kappa [Fix] : 1.000
## lambda [ ML] : 0.751
## lower bound : 0.000, p = 0.022585
## upper bound : 1.000, p = 0.055727
## 95.0% CI : (0.103, NA)
## delta [Fix] : 1.000
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.33659421 0.13968505 2.4097 0.02627 *
## Mass 0.00063677 0.00066595 0.9562 0.35099
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.02915 on 19 degrees of freedom
## Multiple R-squared: 0.04591, Adjusted R-squared: -0.004304
## F-statistic: 0.9143 on 1 and 19 DF, p-value: 0.351pgls_WinFreq_Mass2_3_nosp <- pgls(formula = WinFreq ~ Mass2_3, data = both_comparative_nosp, lambda = "ML")
summary(pgls_WinFreq_Mass2_3_nosp)##
## Call:
## pgls(formula = WinFreq ~ Mass2_3, data = both_comparative_nosp,
## lambda = "ML")
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.061849 -0.008427 0.005194 0.019045 0.054724
##
## Branch length transformations:
##
## kappa [Fix] : 1.000
## lambda [ ML] : 0.746
## lower bound : 0.000, p = 0.023835
## upper bound : 1.000, p = 0.055536
## 95.0% CI : (0.095, NA)
## delta [Fix] : 1.000
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.3137092 0.1528525 2.0524 0.05417 .
## Mass2_3 0.0045054 0.0046532 0.9682 0.34510
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.02906 on 19 degrees of freedom
## Multiple R-squared: 0.04702, Adjusted R-squared: -0.003137
## F-statistic: 0.9375 on 1 and 19 DF, p-value: 0.3451# Check the distribution of the residuals:
par(mfrow=c(2,2))
plot(pgls_CBI_Mass)
par(mfrow=c(2,2))
plot(pgls_CBI_Mass2_3)
par(mfrow=c(2,2))
plot(pgls_CBI_Mass_nosp)
par(mfrow=c(2,2))
plot(pgls_CBI_Mass2_3_nosp)
par(mfrow=c(2,2))
plot(pgls_WinFreq_Mass)
par(mfrow=c(2,2))
plot(pgls_WinFreq_Mass2_3)
par(mfrow=c(2,2))
plot(pgls_WinFreq_Mass_nosp)
par(mfrow=c(2,2))
plot(pgls_WinFreq_Mass2_3_nosp)
This worked great: the \(\lambda\) and p values for the winning frequency are the same as before, and while \(\lambda\) is still running up against the upper bound of 1 for CBI, the fact that it cannot exceed it results in much more sensible p values. We can now proceed to visualizing the results:
# Generate a new Figure 1:
png("/Users/David/Documents/College/2017–18/Winter 2018/FBQ paper/Behav Ecol revisions/Figure_1.png", width = 3600, height = 3600, pointsize = 96)
# setEPS()
# postscript("/Users/David/Documents/College/2017–18/Winter 2018/FBQ_git/Figure_1.eps")
spec_dec <- function(x, k) trimws(format(round(x, k), nsmall = k)) # credit: Jeromy Anglim, Stack Overflow
layout(matrix(c(seq(1,4,1)), nrow=2, ncol=2, byrow=TRUE))
par(mar=c(5, 4, 2, 2))
plot(both_phylo_test$Mass, both_phylo_test$CBI, xlab = "Body mass (g)", ylab = "CBI", main = "Isometric scaling", pch = 16, col = ifelse(both_phylo_test$Mass == both_phylo_test$Mass[both_phylo_test$Species == "Spotted Porcupinefish"], "gray70", "black"))
title("(a)", line = -1, adj = 0, outer = TRUE)
title("(b)", line = -1, outer = TRUE)
title("(c)", line = -20, adj = 0, outer = TRUE)
title("(d)", line = -20, outer = TRUE)
abline(pgls_CBI_Mass, lwd = 5)
abline(pgls_CBI_Mass_nosp, lwd = 5, col = "gray70")
lambda <- pgls_CBI_Mass$param[2]
pseudor2 <- summary(pgls_CBI_Mass)$r.squared
pvalue <- summary(pgls_CBI_Mass)$coef[2,4]
lambda_nosp <- pgls_CBI_Mass_nosp$param[2]
pseudor2_nosp <- summary(pgls_CBI_Mass_nosp)$r.squared
pvalue_nosp <- summary(pgls_CBI_Mass_nosp)$coef[2,4]
text(600, 2.95, cex = 0.8, bquote(lambda == .(spec_dec(lambda, 3))), pos = 4)
text(600, 2.75, cex = 0.8, bquote(italic(R)^2 == .(spec_dec(pseudor2, 3))), pos = 4)
text(600, 2.5, cex = 0.8, bquote(italic(p) == .(spec_dec(pvalue, 3))), pos = 4)
text(600, 1.3, cex = 0.8, bquote(lambda == .(spec_dec(lambda_nosp, 3))), col = "gray70", pos = 4)
text(600, 1.1, cex = 0.8, bquote(italic(R)^2 == .(spec_dec(pseudor2_nosp, 3))), col = "gray70", pos = 4)
text(600, 0.85, cex = 0.8, bquote(italic(p) == .(spec_dec(pvalue_nosp, 3))), col = "gray70", pos = 4)
par(mar=c(5, 4, 2, 2))
plot(both_phylo_test$Mass2_3, both_phylo_test$CBI, xlab = "Body mass raised to 2/3", ylab = "", main = "Allometric scaling", pch = 16, col = ifelse(both_phylo_test$Mass == both_phylo_test$Mass[both_phylo_test$Species == "Spotted Porcupinefish"], "gray70", "black"))
abline(pgls_CBI_Mass2_3, lwd = 5)
abline(pgls_CBI_Mass2_3_nosp, lwd = 5, col = "gray70")
lambda <- pgls_CBI_Mass2_3$param[2]
pseudor2 <- summary(pgls_CBI_Mass2_3)$r.squared
pvalue <- summary(pgls_CBI_Mass2_3)$coef[2,4]
lambda_nosp <- pgls_CBI_Mass2_3_nosp$param[2]
pseudor2_nosp <- summary(pgls_CBI_Mass2_3_nosp)$r.squared
pvalue_nosp <- summary(pgls_CBI_Mass2_3_nosp)$coef[2,4]
text(58, 2.95, cex = 0.8, bquote(lambda == .(spec_dec(lambda, 3))), pos = 4)
text(58, 2.75, cex = 0.8, bquote(italic(R)^2 == .(spec_dec(pseudor2, 3))), pos = 4)
text(58, 2.5, cex = 0.8, bquote(italic(p) == .(spec_dec(pvalue, 3))), pos = 4)
text(58, 2, cex = 0.8, bquote(lambda == .(spec_dec(lambda_nosp, 3))), col = "gray70", pos = 4)
text(58, 1.8, cex = 0.8, bquote(italic(R)^2 == .(spec_dec(pseudor2_nosp, 3))), col = "gray70", pos = 4)
text(58, 1.55, cex = 0.8, bquote(italic(p) == .(spec_dec(pvalue_nosp, 3))), col = "gray70", pos = 4)
par(mar=c(5, 4, 2, 2))
plot(both_phylo_test$Mass, both_phylo_test$WinFreq, xlab = "Body mass (g)", ylab = "Frequency of winning", pch = 16, col = ifelse(both_phylo_test$Mass == both_phylo_test$Mass[both_phylo_test$Species == "Spotted Porcupinefish"], "gray70", "black"))
abline(pgls_WinFreq_Mass, lwd = 5)
abline(pgls_WinFreq_Mass_nosp, lwd = 5, col = "gray70")
lambda <- pgls_WinFreq_Mass$param[2]
pseudor2 <- summary(pgls_WinFreq_Mass)$r.squared
pvalue <- summary(pgls_WinFreq_Mass)$coef[2,4]
lambda_nosp <- pgls_WinFreq_Mass_nosp$param[2]
pseudor2_nosp <- summary(pgls_WinFreq_Mass_nosp)$r.squared
pvalue_nosp <- summary(pgls_WinFreq_Mass_nosp)$coef[2,4]
text(600, 0.7, cex = 0.8, bquote(lambda == .(spec_dec(lambda, 3))), pos = 4)
text(600, 0.7 - 1.2/18, cex = 0.8, bquote(italic(R)^2 == .(spec_dec(pseudor2, 3))), pos = 4)
text(600, 0.7 - 1.2/8, cex = 0.8, bquote(italic(p) == .(spec_dec(pvalue, 3))), pos = 4)
text(600, 0.225, cex = 0.8, bquote(lambda == .(spec_dec(lambda_nosp, 3))), col = "gray70", pos = 4)
text(600, 0.225 - 1.2/18, cex = 0.8, bquote(italic(R)^2 == .(spec_dec(pseudor2_nosp, 3))), col = "gray70", pos = 4)
text(600, 0.225 - 1.2/8, cex = 0.8, bquote(italic(p) == .(spec_dec(pvalue_nosp, 3))), col = "gray70", pos = 4)
par(mar=c(5, 4, 2, 2))
plot(both_phylo_test$Mass2_3, both_phylo_test$WinFreq, xlab = "Body mass raised to 2/3", ylab = "", pch = 16, col = ifelse(both_phylo_test$Mass == both_phylo_test$Mass[both_phylo_test$Species == "Spotted Porcupinefish"], "gray70", "black"))
abline(pgls_WinFreq_Mass2_3, lwd = 5)
abline(pgls_WinFreq_Mass2_3_nosp, lwd = 5, col = "gray70")
lambda <- pgls_WinFreq_Mass2_3$param[2]
pseudor2 <- summary(pgls_WinFreq_Mass2_3)$r.squared
pvalue <- summary(pgls_WinFreq_Mass2_3)$coef[2,4]
lambda_nosp <- pgls_WinFreq_Mass2_3_nosp$param[2]
pseudor2_nosp <- summary(pgls_WinFreq_Mass2_3_nosp)$r.squared
pvalue_nosp <- summary(pgls_WinFreq_Mass2_3_nosp)$coef[2,4]
text(58, 0.98, cex = 0.8, bquote(lambda == .(spec_dec(lambda, 3))), pos = 4)
text(58, 0.98 - 1.2/18, cex = 0.8, bquote(italic(R)^2 == .(spec_dec(pseudor2, 3))), pos = 4)
text(58, 0.98 - 1.2/8, cex = 0.8, bquote(italic(p) == .(spec_dec(pvalue, 3))), pos = 4)
text(58, 0.225, cex = 0.8, bquote(lambda == .(spec_dec(lambda_nosp, 3))), col = "gray70", pos = 4)
text(58, 0.225 - 1.2/18, cex = 0.8, bquote(italic(R)^2 == .(spec_dec(pseudor2_nosp, 3))), col = "gray70", pos = 4)
text(58, 0.225 - 1.2/8, cex = 0.8, bquote(italic(p) == .(spec_dec(pvalue_nosp, 3))), col = "gray70", pos = 4)
dev.off()Group size-only and body size-only models
Below, we conduct another test in which the \(\alpha\) coefficients (individual fighting abilities) are all set to 1, and therefore effectively disregarded. We also evaluate the fit of another model, which includes body size as the only parameter. Note that neither model requires us to repeat the test for body mass raised to 2/3. The group size-only model does not include this parameter at all, while the body size-only model will always yield the same predictions under both mass scalings, since \(m^{\frac{2}{3}}\) is a monotonically increasing function of \(m\).
# Add a new column to both_group with predictions based on group size alone:
for(i in 1:nrow(both_group)) {
if (both_group$Group.A[i] > both_group$Group.B[i]) {
both_group[i,14] <- both_group$Species.A[i]
} else {
both_group[i,14] <- both_group$Species.B[i]
}
}
# Add a new column to both_group with predictions based on body size alone:
for(i in 1:nrow(both_group)) {
if (both_group$BodySize.A[i] > both_group$BodySize.B[i]) {
both_group[i,15] <- both_group$Species.A[i]
} else {
both_group[i,15] <- both_group$Species.B[i]
}
}
colnames(both_group)[14:15] <- c("Group.Size.Only.Prediction","Body.Size.Only.Prediction")
# Test if the group size-only model works:
for(i in 1:nrow(both_group)) {
if (identical(as.character(both_group$Winner[i]), as.character(both_group$Group.Size.Only.Prediction[i])) == TRUE) {
both_group[i,16] <- 1
} else {
both_group[i,16] <- 0
}
}
# Test if the body size-only model works:
for(i in 1:nrow(both_group)) {
if (identical(as.character(both_group$Winner[i]), as.character(both_group$Body.Size.Only.Prediction[i])) == TRUE) {
both_group[i,17] <- 1
} else {
both_group[i,17] <- 0
}
}
colnames(both_group)[16:17] <- c("Group.Size.Only.True", "Body.Size.Only.True")
# Proportion of successful predictions:
sum(both_group$Group.Size.Only.True)/nrow(both_group)## [1] 0.2289157sum(both_group$Body.Size.Only.True)/nrow(both_group)## [1] 0.6987952# Test significance:
dbinom(sum(both_group$Group.Size.Only.True), size = nrow(both_group), prob = 1/2)## [1] 2.643039e-07dbinom(sum(both_group$Body.Size.Only.True), size = nrow(both_group), prob = 1/2)## [1] 0.0001118917In what cases did the body size-only model work but the next best model (Lanchester’s linear law) did not?
both_group[both_group$Body.Size.Only.True == 1 & both_group$Linear.True == 0,]## Species.A Species.B Group.A Group.B
## 59 Damselfish Scissortail Sergeant 1 2
## 87 Damselfish Sixbar Wrasse 1 2
## 246 Damselfish Juvenile Parrotfish 1 20
## 259 Regal Angelfish Brown Tang 1 2
## 273 Orange-lined Triggerfish Threeband Pennantfish 1 2
## 276 Orange-lined Triggerfish Brown Tang 1 2
## 324 Damselfish Scissortail Sergeant 1 3
## 329 Damselfish Sixbar Wrasse 1 2
## 487 Damselfish Juvenile Parrotfish 1 7
## 515 Damselfish Juvenile Parrotfish 1 14
## 554 Damselfish Scissortail Sergeant 1 2
## Winner BodySize.A BodySize.B Linear.Law.Prediction
## 59 Damselfish 38.88513 25.35153 Scissortail Sergeant
## 87 Damselfish 38.88513 22.72689 Sixbar Wrasse
## 246 Damselfish 38.88513 7.49326 Juvenile Parrotfish
## 259 Regal Angelfish 128.17052 78.11946 Brown Tang
## 273 Orange-lined Triggerfish 108.96596 87.65761 Threeband Pennantfish
## 276 Orange-lined Triggerfish 108.96596 78.11946 Brown Tang
## 324 Damselfish 38.88513 25.35153 Scissortail Sergeant
## 329 Damselfish 38.88513 22.72689 Sixbar Wrasse
## 487 Damselfish 38.88513 7.49326 Juvenile Parrotfish
## 515 Damselfish 38.88513 7.49326 Juvenile Parrotfish
## 554 Damselfish 38.88513 25.35153 Scissortail Sergeant
## Square.Law.Prediction Linear.True Square.True
## 59 Scissortail Sergeant 0 0
## 87 Sixbar Wrasse 0 0
## 246 Juvenile Parrotfish 0 0
## 259 Brown Tang 0 0
## 273 Threeband Pennantfish 0 0
## 276 Brown Tang 0 0
## 324 Scissortail Sergeant 0 0
## 329 Sixbar Wrasse 0 0
## 487 Juvenile Parrotfish 0 0
## 515 Juvenile Parrotfish 0 0
## 554 Scissortail Sergeant 0 0
## Diet.A Diet.B
## 59 Omnivore Omnivore
## 87 Omnivore Benthic Invertebrate Consumer
## 246 Omnivore Herbivore/Detritivore
## 259 Benthic Invertebrate Consumer Herbivore/Detritivore
## 273 Omnivore Corallivore
## 276 Omnivore Herbivore/Detritivore
## 324 Omnivore Omnivore
## 329 Omnivore Benthic Invertebrate Consumer
## 487 Omnivore Herbivore/Detritivore
## 515 Omnivore Herbivore/Detritivore
## 554 Omnivore Omnivore
## Group.Size.Only.Prediction Body.Size.Only.Prediction
## 59 Scissortail Sergeant Damselfish
## 87 Sixbar Wrasse Damselfish
## 246 Juvenile Parrotfish Damselfish
## 259 Brown Tang Regal Angelfish
## 273 Threeband Pennantfish Orange-lined Triggerfish
## 276 Brown Tang Orange-lined Triggerfish
## 324 Scissortail Sergeant Damselfish
## 329 Sixbar Wrasse Damselfish
## 487 Juvenile Parrotfish Damselfish
## 515 Juvenile Parrotfish Damselfish
## 554 Scissortail Sergeant Damselfish
## Group.Size.Only.True Body.Size.Only.True
## 59 0 1
## 87 0 1
## 246 0 1
## 259 0 1
## 273 0 1
## 276 0 1
## 324 0 1
## 329 0 1
## 487 0 1
## 515 0 1
## 554 0 18 out of the 11 interactions in this category involved damselfish. What happens if we remove them from the dataset?
both_group3 <- both_group[both_group$Species.A != "Damselfish" & both_group$Species.B != "Damselfish", ]
nrow(both_group3)## [1] 43# Proportion of successful predictions:
sum(both_group3$Linear.True)/nrow(both_group3)## [1] 0.8372093sum(both_group3$Square.True)/nrow(both_group3)## [1] 0.7209302sum(both_group3$Group.Size.Only.True)/nrow(both_group3)## [1] 0.255814sum(both_group3$Body.Size.Only.True)/nrow(both_group3)## [1] 0.7906977# Test significance:
dbinom(sum(both_group3$Linear.True), size = nrow(both_group3), prob = 1/2)## [1] 3.663458e-06dbinom(sum(both_group3$Square.True), size = nrow(both_group3), prob = 1/2)## [1] 0.001743806dbinom(sum(both_group3$Group.Size.Only.True), size = nrow(both_group3), prob = 1/2)## [1] 0.0006539272dbinom(sum(both_group3$Body.Size.Only.True), size = nrow(both_group3), prob = 1/2)## [1] 6.411051e-05# With mass raised to the two-thirds power. Remember, this can change the results for the
# linear and square laws, but not for the body size-only or group size-only models.
both_group4 <- both_group2[both_group2$Species.A != "Damselfish" & both_group2$Species.B != "Damselfish" ,]
nrow(both_group4)## [1] 43# Proportion of successful predictions:
sum(both_group4$Linear.True)/nrow(both_group4)## [1] 0.7906977sum(both_group4$Square.True)/nrow(both_group4)## [1] 0.5348837# Test significance:
dbinom(sum(both_group4$Linear.True), size = nrow(both_group4), prob = 1/2)## [1] 6.411051e-05dbinom(sum(both_group4$Square.True), size = nrow(both_group4), prob = 1/2)## [1] 0.1092038In the next step, we calculate the statistical power of the linear law, square law, and body size-only model tests. All of the below is based on http://www.calvin.edu/~scofield/courses/m343/F15/handouts/binomialTestPower.pdf
# First, find the rejection region for n = 83 (the number of data points in both_group)
# under the null assumption (probability of successful prediction = 0.5) and significance
# level alpha = 0.05:
qbinom(0.025, nrow(both_group), 0.5) # 33## [1] 33qbinom(0.975, nrow(both_group), 0.5) # 50## [1] 50# Since P(33 \leq x \geq 50) is probably not equal to exactly 0.95, we need to find out
# if the sum of the values below is larger or smaller than this significance threshold:
pbinom(33, nrow(both_group), 0.5) + (1 - pbinom(50, nrow(both_group), 0.5))## [1] 0.06297226pbinom(32, nrow(both_group), 0.5) + (1 - pbinom(51, nrow(both_group), 0.5))## [1] 0.03752953# Let's use 32 and 51 as the boundaries of the rejection region.
# Compute power as 1 - beta:
1 - (pbinom(51, nrow(both_group), sum(both_group$Body.Size.Only.True)/nrow(both_group)) -
pbinom(32, nrow(both_group), sum(both_group$Body.Size.Only.True)/nrow(both_group)))## [1] 0.9377759Compute the maximum likelihood estimate under the binomial and use it to compute AIC scores:
# -2*log(L^hat) for the body size-only model:
-2*log(dbinom(sum(both_group$Body.Size.Only.True),
nrow(both_group),
sum(both_group$Body.Size.Only.True)/nrow(both_group)))## [1] 4.705887# -2*log(L^hat) for the linear law:
-2*log(dbinom(sum(both_group$Linear.True),
nrow(both_group),
sum(both_group$Linear.True)/nrow(both_group)))## [1] 4.765549Finally, see how the body size-only and group-size only models perform when fitted to the datasets subsampled by diet:
test_omni## Species.A Species.B Group.A Group.B
## 59 Damselfish Scissortail Sergeant 1 2
## 113 Moorish Idol Lemonpeel Angelfish 1 2
## 132 Damselfish Lemonpeel Angelfish 1 2
## 187 Dusky Angelfish Scissortail Sergeant 1 3
## 288 Orange-lined Triggerfish Lemonpeel Angelfish 1 2
## 324 Damselfish Scissortail Sergeant 1 3
## 353 Orange-lined Triggerfish Lemonpeel Angelfish 1 3
## 355 Damselfish Lemonpeel Angelfish 1 3
## 357 Damselfish Lemonpeel Angelfish 1 3
## 358 Damselfish Lemonpeel Angelfish 2 1
## 363 Orange-lined Triggerfish Lemonpeel Angelfish 1 3
## 366 Moorish Idol Lemonpeel Angelfish 1 3
## 368 Damselfish Lemonpeel Angelfish 1 3
## 373 Damselfish Lemonpeel Angelfish 1 3
## 440 Damselfish Lemonpeel Angelfish 1 2
## 554 Damselfish Scissortail Sergeant 1 2
## Winner BodySize.A BodySize.B
## 59 Damselfish 38.88513 25.351526
## 113 Moorish Idol 171.38610 7.490788
## 132 Damselfish 38.88513 7.490788
## 187 Scissortail Sergeant 10.99104 25.351526
## 288 Orange-lined Triggerfish 108.96596 7.490788
## 324 Damselfish 38.88513 25.351526
## 353 Orange-lined Triggerfish 108.96596 7.490788
## 355 Damselfish 38.88513 7.490788
## 357 Damselfish 38.88513 7.490788
## 358 Damselfish 38.88513 7.490788
## 363 Orange-lined Triggerfish 108.96596 7.490788
## 366 Moorish Idol 171.38610 7.490788
## 368 Damselfish 38.88513 7.490788
## 373 Damselfish 38.88513 7.490788
## 440 Damselfish 38.88513 7.490788
## 554 Damselfish 38.88513 25.351526
## Linear.Law.Prediction Square.Law.Prediction Linear.True
## 59 Scissortail Sergeant Scissortail Sergeant 0
## 113 Moorish Idol Moorish Idol 1
## 132 Damselfish Damselfish 1
## 187 Scissortail Sergeant Scissortail Sergeant 1
## 288 Orange-lined Triggerfish Orange-lined Triggerfish 1
## 324 Scissortail Sergeant Scissortail Sergeant 0
## 353 Orange-lined Triggerfish Orange-lined Triggerfish 1
## 355 Damselfish Lemonpeel Angelfish 1
## 357 Damselfish Lemonpeel Angelfish 1
## 358 Damselfish Damselfish 1
## 363 Orange-lined Triggerfish Orange-lined Triggerfish 1
## 366 Moorish Idol Moorish Idol 1
## 368 Damselfish Lemonpeel Angelfish 1
## 373 Damselfish Lemonpeel Angelfish 1
## 440 Damselfish Damselfish 1
## 554 Scissortail Sergeant Scissortail Sergeant 0
## Square.True Diet.A Diet.B
## 59 0 Omnivore Omnivore
## 113 1 Omnivore Omnivore
## 132 1 Omnivore Omnivore
## 187 1 Omnivore Omnivore
## 288 1 Omnivore Omnivore
## 324 0 Omnivore Omnivore
## 353 1 Omnivore Omnivore
## 355 0 Omnivore Omnivore
## 357 0 Omnivore Omnivore
## 358 1 Omnivore Omnivore
## 363 1 Omnivore Omnivore
## 366 1 Omnivore Omnivore
## 368 0 Omnivore Omnivore
## 373 0 Omnivore Omnivore
## 440 1 Omnivore Omnivore
## 554 0 Omnivore Omnivoretest_omni2## Species.A Species.B Group.A Group.B
## 59 Damselfish Scissortail Sergeant 1 2
## 113 Moorish Idol Lemonpeel Angelfish 1 2
## 132 Damselfish Lemonpeel Angelfish 1 2
## 187 Dusky Angelfish Scissortail Sergeant 1 3
## 288 Orange-lined Triggerfish Lemonpeel Angelfish 1 2
## 324 Damselfish Scissortail Sergeant 1 3
## 353 Orange-lined Triggerfish Lemonpeel Angelfish 1 3
## 355 Damselfish Lemonpeel Angelfish 1 3
## 357 Damselfish Lemonpeel Angelfish 1 3
## 358 Damselfish Lemonpeel Angelfish 2 1
## 363 Orange-lined Triggerfish Lemonpeel Angelfish 1 3
## 366 Moorish Idol Lemonpeel Angelfish 1 3
## 368 Damselfish Lemonpeel Angelfish 1 3
## 373 Damselfish Lemonpeel Angelfish 1 3
## 440 Damselfish Lemonpeel Angelfish 1 2
## 554 Damselfish Scissortail Sergeant 1 2
## Winner BodySize.A BodySize.B
## 59 Damselfish 11.477723 8.629840
## 113 Moorish Idol 30.854397 3.828409
## 132 Damselfish 11.477723 3.828409
## 187 Scissortail Sergeant 4.943401 8.629840
## 288 Orange-lined Triggerfish 22.813604 3.828409
## 324 Damselfish 11.477723 8.629840
## 353 Orange-lined Triggerfish 22.813604 3.828409
## 355 Damselfish 11.477723 3.828409
## 357 Damselfish 11.477723 3.828409
## 358 Damselfish 11.477723 3.828409
## 363 Orange-lined Triggerfish 22.813604 3.828409
## 366 Moorish Idol 30.854397 3.828409
## 368 Damselfish 11.477723 3.828409
## 373 Damselfish 11.477723 3.828409
## 440 Damselfish 11.477723 3.828409
## 554 Damselfish 11.477723 8.629840
## Linear.Law.Prediction Square.Law.Prediction Linear.True
## 59 Scissortail Sergeant Scissortail Sergeant 0
## 113 Moorish Idol Moorish Idol 1
## 132 Damselfish Lemonpeel Angelfish 1
## 187 Scissortail Sergeant Scissortail Sergeant 1
## 288 Orange-lined Triggerfish Orange-lined Triggerfish 1
## 324 Scissortail Sergeant Scissortail Sergeant 0
## 353 Orange-lined Triggerfish Lemonpeel Angelfish 1
## 355 Lemonpeel Angelfish Lemonpeel Angelfish 0
## 357 Lemonpeel Angelfish Lemonpeel Angelfish 0
## 358 Damselfish Damselfish 1
## 363 Orange-lined Triggerfish Lemonpeel Angelfish 1
## 366 Moorish Idol Lemonpeel Angelfish 1
## 368 Lemonpeel Angelfish Lemonpeel Angelfish 0
## 373 Lemonpeel Angelfish Lemonpeel Angelfish 0
## 440 Damselfish Lemonpeel Angelfish 1
## 554 Scissortail Sergeant Scissortail Sergeant 0
## Square.True Diet.A Diet.B
## 59 0 Omnivore Omnivore
## 113 1 Omnivore Omnivore
## 132 0 Omnivore Omnivore
## 187 1 Omnivore Omnivore
## 288 1 Omnivore Omnivore
## 324 0 Omnivore Omnivore
## 353 0 Omnivore Omnivore
## 355 0 Omnivore Omnivore
## 357 0 Omnivore Omnivore
## 358 1 Omnivore Omnivore
## 363 0 Omnivore Omnivore
## 366 0 Omnivore Omnivore
## 368 0 Omnivore Omnivore
## 373 0 Omnivore Omnivore
## 440 0 Omnivore Omnivore
## 554 0 Omnivore Omnivoretest_herb## Species.A Species.B Group.A Group.B
## 116 Brown Tang Juvenile Parrotfish 1 4
## 182 Dark-capped Parrotfish Juvenile Parrotfish 1 10
## 264 Bullethead Parrotfish Brown Tang 1 2
## 282 Dark-capped Parrotfish Juvenile Parrotfish 1 20
## Winner BodySize.A BodySize.B Linear.Law.Prediction
## 116 Brown Tang 78.11946 7.49326 Brown Tang
## 182 Dark-capped Parrotfish 311.91661 7.49326 Dark-capped Parrotfish
## 264 Bullethead Parrotfish 339.38694 78.11946 Bullethead Parrotfish
## 282 Dark-capped Parrotfish 311.91661 7.49326 Dark-capped Parrotfish
## Square.Law.Prediction Linear.True Square.True Diet.A
## 116 Juvenile Parrotfish 1 0 Herbivore/Detritivore
## 182 Juvenile Parrotfish 1 0 Herbivore/Detritivore
## 264 Bullethead Parrotfish 1 1 Herbivore/Detritivore
## 282 Juvenile Parrotfish 1 0 Herbivore/Detritivore
## Diet.B
## 116 Herbivore/Detritivore
## 182 Herbivore/Detritivore
## 264 Herbivore/Detritivore
## 282 Herbivore/Detritivoretest_herb2## Species.A Species.B Group.A Group.B
## 116 Brown Tang Juvenile Parrotfish 1 4
## 182 Dark-capped Parrotfish Juvenile Parrotfish 1 10
## 264 Bullethead Parrotfish Brown Tang 1 2
## 282 Dark-capped Parrotfish Juvenile Parrotfish 1 20
## Winner BodySize.A BodySize.B Linear.Law.Prediction
## 116 Brown Tang 18.27425 3.829251 Brown Tang
## 182 Dark-capped Parrotfish 45.99306 3.829251 Dark-capped Parrotfish
## 264 Bullethead Parrotfish 48.65529 18.274246 Bullethead Parrotfish
## 282 Dark-capped Parrotfish 45.99306 3.829251 Juvenile Parrotfish
## Square.Law.Prediction Linear.True Square.True Diet.A
## 116 Juvenile Parrotfish 1 0 Herbivore/Detritivore
## 182 Juvenile Parrotfish 1 0 Herbivore/Detritivore
## 264 Brown Tang 1 0 Herbivore/Detritivore
## 282 Juvenile Parrotfish 0 0 Herbivore/Detritivore
## Diet.B
## 116 Herbivore/Detritivore
## 182 Herbivore/Detritivore
## 264 Herbivore/Detritivore
## 282 Herbivore/Detritivore# The function below adds columns "Group.Size.Only.Prediction",
# "Body.Size.Only.Prediction", "Group.Size.Only.True", and "Body.Size.Only.True" to an
# arbitrary data frame containing columns named "Group.A", "Group.B", "BodySizeA",
# "BodySizeB", and "Winner":
test_extra_two_models <- function(dataframe) {
n <- ncol(dataframe)
for(i in 1:nrow(dataframe)) {
if (dataframe$Group.A[i] > dataframe$Group.B[i]) {
dataframe[i, n+1] <- dataframe$Species.A[i]
} else {
dataframe[i, n+1] <- dataframe$Species.B[i]
}
if (dataframe$BodySize.A[i] > dataframe$BodySize.B[i]) {
dataframe[i, n+2] <- dataframe$Species.A[i]
} else {
dataframe[i, n+2] <- dataframe$Species.B[i]
}
if (identical(as.character(dataframe$Winner[i]), as.character(dataframe[i, n+1])) == TRUE) {
dataframe[i, n+3] <- 1
} else {
dataframe[i, n+3] <- 0
}
if (identical(as.character(dataframe$Winner[i]), as.character(dataframe[i, n+2])) == TRUE) {
dataframe[i, n+4] <- 1
} else {
dataframe[i, n+4] <- 0
}
}
colnames(dataframe)[(n+1):(n+4)] <- c("Group.Size.Only.Prediction",
"Body.Size.Only.Prediction",
"Group.Size.Only.True",
"Body.Size.Only.True")
return(dataframe)
}
test_omni <- test_extra_two_models(test_omni)
test_herb <- test_extra_two_models(test_herb)
# Proportion of successful predictions:
sum(test_omni$Group.Size.Only.True)/nrow(test_omni)## [1] 0.125sum(test_omni$Body.Size.Only.True)/nrow(test_omni)## [1] 1sum(test_herb$Group.Size.Only.True)/nrow(test_herb)## [1] 0sum(test_herb$Body.Size.Only.True)/nrow(test_herb)## [1] 1# Test significance:
dbinom(sum(test_omni$Group.Size.Only.True), size = nrow(test_omni), prob = 1/2)## [1] 0.001831055dbinom(sum(test_omni$Body.Size.Only.True), size = nrow(test_omni), prob = 1/2)## [1] 1.525879e-05dbinom(sum(test_herb$Group.Size.Only.True), size = nrow(test_herb), prob = 1/2)## [1] 0.0625dbinom(sum(test_herb$Body.Size.Only.True), size = nrow(test_herb), prob = 1/2)## [1] 0.0625Live/video differences
We have already alculated the properties of matrices based on live-only and video-only observations. In the next step, we compare the hierarchies:
# Subsample each taxon-by-taxon dominance matrix down to the shared set of 15 species
nojp_live_shared <- exclude(nojp_live, setdiff(rownames(nojp_live), rownames(nojp_scored)))
nojp_scored_shared <- exclude(nojp_scored, setdiff(rownames(nojp_scored), rownames(nojp_live)))
nojp_both_shared <- exclude(nojp_both, rownames(nojp_both)[!rownames(nojp_both) %in% intersect(rownames(nojp_live), rownames(nojp_scored))])
# Recalculate the CBI for the 15 species based on the live-only and video-only datasets:
live_subsampled_cbi <- cbi_hierarchy_new(nojp_live_shared)
colnames(live_subsampled_cbi)[2] <- "Live"
video_subsampled_cbi <- cbi_hierarchy_new(nojp_scored_shared)
colnames(video_subsampled_cbi)[2] <- "Video"
subsampled_cbi <- merge(live_subsampled_cbi, video_subsampled_cbi, by = "Species")
both_subsampled_cbi <- cbi_hierarchy_new(nojp_both_shared)
subsampled_cbi <- merge(subsampled_cbi, both_subsampled_cbi, by = "Species")
colnames(subsampled_cbi)[4] <- "Both"
# Recalculate the frequency of winning for the 15 species based on the live-only and
# video-only datasets:
live_subsampled_winfreq <- win_freq(sheet_ooo5)[win_freq(sheet_ooo5)$Species %in% intersect(win_freq(sheet_ooo5)$Species, win_freq(scored_ooo5)$Species),]
live_subsampled_winfreq <- live_subsampled_winfreq[order(live_subsampled_winfreq$Species),]
colnames(live_subsampled_winfreq)[2] <- "Live"
video_subsampled_winfreq <- win_freq(scored_ooo5)[win_freq(scored_ooo5)$Species %in% intersect(win_freq(sheet_ooo5)$Species, win_freq(scored_ooo5)$Species),]
video_subsampled_winfreq <- video_subsampled_winfreq[order(video_subsampled_winfreq$Species),]
colnames(video_subsampled_winfreq)[2] <- "Video"
subsampled_winfreq <- merge(live_subsampled_winfreq, video_subsampled_winfreq, by = "Species")
subsampled_winfreq <- merge(subsampled_winfreq, win_freq(both_ooo5), by = "Species")
colnames(subsampled_winfreq)[4] <- "Both"
if(!require("RcmdrMisc")) {install.packages("RcmdrMisc")} # Beware: has many dependencies!## Loading required package: RcmdrMisc## Loading required package: car## Loading required package: carData## Loading required package: sandwichlibrary(RcmdrMisc)
rcorr.adjust(as.matrix(subsampled_cbi[,-1]), type = "spearman")##
## Spearman correlations:
## Live Video Both
## Live 1.0000 0.2898 0.9357
## Video 0.2898 1.0000 0.3792
## Both 0.9357 0.3792 1.0000
##
## Number of observations: 15
##
## Pairwise two-sided p-values:
## Live Video Both
## Live 0.2948 <.0001
## Video 0.2948 0.1633
## Both <.0001 0.1633
##
## Adjusted p-values (Holm's method)
## Live Video Both
## Live 0.3266 <.0001
## Video 0.3266 0.3266
## Both <.0001 0.3266rcorr.adjust(as.matrix(subsampled_cbi[,-1]), type = "pearson")##
## Pearson correlations:
## Live Video Both
## Live 1.0000 0.6703 0.8643
## Video 0.6703 1.0000 0.8914
## Both 0.8643 0.8914 1.0000
##
## Number of observations: 15
##
## Pairwise two-sided p-values:
## Live Video Both
## Live 0.0062 <.0001
## Video 0.0062 <.0001
## Both <.0001 <.0001
##
## Adjusted p-values (Holm's method)
## Live Video Both
## Live 0.0062 <.0001
## Video 0.0062 <.0001
## Both <.0001 <.0001rcorr.adjust(as.matrix(subsampled_winfreq[,-1]), type = "spearman")##
## Spearman correlations:
## Live Video Both
## Live 1.0000 0.6188 0.8959
## Video 0.6188 1.0000 0.8307
## Both 0.8959 0.8307 1.0000
##
## Number of observations: 15
##
## Pairwise two-sided p-values:
## Live Video Both
## Live 0.0139 <.0001
## Video 0.0139 0.0001
## Both <.0001 0.0001
##
## Adjusted p-values (Holm's method)
## Live Video Both
## Live 0.0139 <.0001
## Video 0.0139 0.0003
## Both <.0001 0.0003rcorr.adjust(as.matrix(subsampled_winfreq[,-1]), type = "pearson")##
## Pearson correlations:
## Live Video Both
## Live 1.0000 0.6208 0.9133
## Video 0.6208 1.0000 0.8731
## Both 0.9133 0.8731 1.0000
##
## Number of observations: 15
##
## Pairwise two-sided p-values:
## Live Video Both
## Live 0.0135 <.0001
## Video 0.0135 <.0001
## Both <.0001 <.0001
##
## Adjusted p-values (Holm's method)
## Live Video Both
## Live 0.0135 <.0001
## Video 0.0135 <.0001
## Both <.0001 <.0001Re-test the relationship between body size and individual fighting ability using live-only and video-only data:
# Appropriately subsample 'fishtree':
fishtree_live <- drop.tip(fishtree, fishtree$tip.label[!fishtree$tip.label %in% rownames(nojp_live)])
fishtree_live_nosp <- drop.tip(fishtree, c("Spotted Porcupinefish", fishtree$tip.label[!fishtree$tip.label %in% rownames(nojp_live)]))
fishtree_video <- drop.tip(fishtree, fishtree$tip.label[!fishtree$tip.label %in% rownames(nojp_scored)])
library(caper)
live_comparative <- comparative.data(data = sheet_phylo_test, phy = fishtree_live, names.col = "Species", vcv = TRUE, warn.dropped = TRUE)
live_comparative_nosp <- comparative.data(data = sheet_phylo_test_nosp, phy = fishtree_live_nosp, names.col = "Species", vcv = TRUE, warn.dropped = TRUE)
video_comparative <- comparative.data(data = scored_phylo_test, phy = fishtree_video, names.col = "Species", vcv = TRUE, warn.dropped = TRUE)
live_pgls_CBI_Mass <- pgls(formula = CBI ~ Mass, data = live_comparative, lambda = "ML")
summary(live_pgls_CBI_Mass)##
## Call:
## pgls(formula = CBI ~ Mass, data = live_comparative, lambda = "ML")
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.12134 -0.09154 -0.02659 0.07927 0.14682
##
## Branch length transformations:
##
## kappa [Fix] : 1.000
## lambda [ ML] : 1.000
## lower bound : 0.000, p = 0.00051509
## upper bound : 1.000, p = 1
## 95.0% CI : (0.803, NA)
## delta [Fix] : 1.000
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.04492043 0.41860472 2.4962 0.02192 *
## Mass -0.00048654 0.00074678 -0.6515 0.52251
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.09167 on 19 degrees of freedom
## Multiple R-squared: 0.02185, Adjusted R-squared: -0.02963
## F-statistic: 0.4245 on 1 and 19 DF, p-value: 0.5225live_pgls_CBI_Mass2_3 <- pgls(formula = CBI ~ Mass2_3, data = live_comparative, lambda = "ML"); summary(live_pgls_CBI_Mass2_3)##
## Call:
## pgls(formula = CBI ~ Mass2_3, data = live_comparative, lambda = "ML")
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.14217 -0.08418 -0.02553 0.07819 0.14423
##
## Branch length transformations:
##
## kappa [Fix] : 1.000
## lambda [ ML] : 1.000
## lower bound : 0.000, p = 0.00083993
## upper bound : 1.000, p = 1
## 95.0% CI : (0.792, NA)
## delta [Fix] : 1.000
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.0018370 0.4523589 2.2147 0.0392 *
## Mass2_3 -0.0018632 0.0078794 -0.2365 0.8156
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.09256 on 19 degrees of freedom
## Multiple R-squared: 0.002934, Adjusted R-squared: -0.04954
## F-statistic: 0.05591 on 1 and 19 DF, p-value: 0.8156live_pgls_WinFreq_Mass <- pgls(formula = WinFreq ~ Mass, data = live_comparative, lambda = "ML"); summary(live_pgls_WinFreq_Mass)##
## Call:
## pgls(formula = WinFreq ~ Mass, data = live_comparative, lambda = "ML")
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.037100 -0.016062 0.006017 0.029263 0.070629
##
## Branch length transformations:
##
## kappa [Fix] : 1.000
## lambda [ ML] : 0.663
## lower bound : 0.000, p = 0.10993
## upper bound : 1.000, p = 0.01489
## 95.0% CI : (NA, 0.973)
## delta [Fix] : 1.000
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.41376043 0.13109752 3.1561 0.005201 **
## Mass 0.00020690 0.00027338 0.7568 0.458435
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.03139 on 19 degrees of freedom
## Multiple R-squared: 0.02926, Adjusted R-squared: -0.02183
## F-statistic: 0.5728 on 1 and 19 DF, p-value: 0.4584live_pgls_WinFreq_Mass2_3 <- pgls(formula = WinFreq ~ Mass2_3, data = live_comparative, lambda = "ML"); summary(live_pgls_WinFreq_Mass2_3)##
## Call:
## pgls(formula = WinFreq ~ Mass2_3, data = live_comparative, lambda = "ML")
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.067679 -0.011879 -0.003656 0.017396 0.048635
##
## Branch length transformations:
##
## kappa [Fix] : 1.000
## lambda [ ML] : 0.631
## lower bound : 0.000, p = 0.13084
## upper bound : 1.000, p = 0.011061
## 95.0% CI : (NA, 0.966)
## delta [Fix] : 1.000
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.3593571 0.1397527 2.5714 0.01869 *
## Mass2_3 0.0033572 0.0028881 1.1624 0.25946
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.03047 on 19 degrees of freedom
## Multiple R-squared: 0.06639, Adjusted R-squared: 0.01726
## F-statistic: 1.351 on 1 and 19 DF, p-value: 0.2595live_pgls_CBI_Mass_nosp <- pgls(formula = CBI ~ Mass, data = live_comparative_nosp, lambda = "ML"); summary(live_pgls_CBI_Mass_nosp)##
## Call:
## pgls(formula = CBI ~ Mass, data = live_comparative_nosp, lambda = "ML")
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.140248 -0.004507 0.026080 0.089520 0.192394
##
## Branch length transformations:
##
## kappa [Fix] : 1.000
## lambda [ ML] : 1.000
## lower bound : 0.000, p = 0.00037668
## upper bound : 1.000, p = 1
## 95.0% CI : (0.787, NA)
## delta [Fix] : 1.000
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.8128104 0.4256586 1.9095 0.07226 .
## Mass 0.0017328 0.0016737 1.0353 0.31424
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.0873 on 18 degrees of freedom
## Multiple R-squared: 0.0562, Adjusted R-squared: 0.003767
## F-statistic: 1.072 on 1 and 18 DF, p-value: 0.3142live_pgls_CBI_Mass2_3_nosp <- pgls(formula = CBI ~ Mass2_3, data = live_comparative_nosp, lambda = "ML"); summary(live_pgls_CBI_Mass2_3_nosp)##
## Call:
## pgls(formula = CBI ~ Mass2_3, data = live_comparative_nosp, lambda = "ML")
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.140775 -0.004609 0.022412 0.089080 0.188876
##
## Branch length transformations:
##
## kappa [Fix] : 1.000
## lambda [ ML] : 1.000
## lower bound : 0.000, p = 0.00039991
## upper bound : 1.000, p = 1
## 95.0% CI : (0.785, NA)
## delta [Fix] : 1.000
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.750749 0.450167 1.6677 0.1127
## Mass2_3 0.012323 0.011456 1.0756 0.2963
##
## Residual standard error: 0.0871 on 18 degrees of freedom
## Multiple R-squared: 0.06039, Adjusted R-squared: 0.008192
## F-statistic: 1.157 on 1 and 18 DF, p-value: 0.2963live_pgls_WinFreq_Mass_nosp <- pgls(formula = WinFreq ~ Mass, data = live_comparative_nosp, lambda = "ML"); summary(live_pgls_WinFreq_Mass_nosp)##
## Call:
## pgls(formula = WinFreq ~ Mass, data = live_comparative_nosp,
## lambda = "ML")
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.014123 0.002685 0.014639 0.036257 0.064441
##
## Branch length transformations:
##
## kappa [Fix] : 1.000
## lambda [ ML] : 0.647
## lower bound : 0.000, p = 0.079667
## upper bound : 1.000, p = 0.0081626
## 95.0% CI : (NA, 0.959)
## delta [Fix] : 1.000
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.31444575 0.14286027 2.2011 0.04102 *
## Mass 0.00116556 0.00070943 1.6430 0.11775
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.03031 on 18 degrees of freedom
## Multiple R-squared: 0.1304, Adjusted R-squared: 0.0821
## F-statistic: 2.699 on 1 and 18 DF, p-value: 0.1177live_pgls_WinFreq_Mass2_3_nosp <- pgls(formula = WinFreq ~ Mass2_3, data = live_comparative_nosp, lambda = "ML"); summary(live_pgls_WinFreq_Mass2_3_nosp)##
## Call:
## pgls(formula = WinFreq ~ Mass2_3, data = live_comparative_nosp,
## lambda = "ML")
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.061269 -0.013488 -0.002605 0.015575 0.042127
##
## Branch length transformations:
##
## kappa [Fix] : 1.000
## lambda [ ML] : 0.643
## lower bound : 0.000, p = 0.079968
## upper bound : 1.000, p = 0.0076034
## 95.0% CI : (NA, 0.957)
## delta [Fix] : 1.000
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.2629693 0.1564604 1.6807 0.11008
## Mass2_3 0.0087316 0.0049639 1.7590 0.09556 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.02997 on 18 degrees of freedom
## Multiple R-squared: 0.1467, Adjusted R-squared: 0.09928
## F-statistic: 3.094 on 1 and 18 DF, p-value: 0.09556video_pgls_CBI_Mass <- pgls(formula = CBI ~ Mass, data = video_comparative, lambda = "ML")
summary(video_pgls_CBI_Mass)##
## Call:
## pgls(formula = CBI ~ Mass, data = video_comparative, lambda = "ML")
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.12487 -0.07136 -0.04121 0.01673 0.10692
##
## Branch length transformations:
##
## kappa [Fix] : 1.000
## lambda [ ML] : 0.049
## lower bound : 0.000, p = 0.87521
## upper bound : 1.000, p = 0.00084113
## 95.0% CI : (NA, 0.719)
## delta [Fix] : 1.000
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.6174704 0.2932005 2.1060 0.05374 .
## Mass -0.0014836 0.0020582 -0.7208 0.48287
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.07338 on 14 degrees of freedom
## Multiple R-squared: 0.03579, Adjusted R-squared: -0.03308
## F-statistic: 0.5196 on 1 and 14 DF, p-value: 0.4829video_pgls_CBI_Mass2_3 <- pgls(formula = CBI ~ Mass2_3, data = video_comparative, lambda = "ML"); summary(video_pgls_CBI_Mass2_3)##
## Call:
## pgls(formula = CBI ~ Mass2_3, data = video_comparative, lambda = "ML")
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.11810 -0.05457 -0.01404 0.03946 0.10942
##
## Branch length transformations:
##
## kappa [Fix] : 1.000
## lambda [ ML] : 0.054
## lower bound : 0.000, p = 0.86644
## upper bound : 1.000, p = 0.00091886
## 95.0% CI : (NA, 0.725)
## delta [Fix] : 1.000
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.6628957 0.3641409 1.8204 0.09013 .
## Mass2_3 -0.0099129 0.0152288 -0.6509 0.52563
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.07362 on 14 degrees of freedom
## Multiple R-squared: 0.02938, Adjusted R-squared: -0.03995
## F-statistic: 0.4237 on 1 and 14 DF, p-value: 0.5256video_pgls_WinFreq_Mass <- pgls(formula = WinFreq ~ Mass, data = video_comparative, lambda = "ML"); summary(video_pgls_WinFreq_Mass)##
## Call:
## pgls(formula = WinFreq ~ Mass, data = video_comparative, lambda = "ML")
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.035535 -0.011008 0.001371 0.023994 0.063956
##
## Branch length transformations:
##
## kappa [Fix] : 1.000
## lambda [ ML] : 0.536
## lower bound : 0.000, p = 0.21789
## upper bound : 1.000, p = 0.15603
## 95.0% CI : (NA, NA)
## delta [Fix] : 1.000
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.28375197 0.13948795 2.0342 0.06133 .
## Mass -0.00043991 0.00083487 -0.5269 0.60650
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.02897 on 14 degrees of freedom
## Multiple R-squared: 0.01945, Adjusted R-squared: -0.05059
## F-statistic: 0.2776 on 1 and 14 DF, p-value: 0.6065video_pgls_WinFreq_Mass2_3 <- pgls(formula = WinFreq ~ Mass2_3, data = video_comparative, lambda = "ML"); summary(video_pgls_WinFreq_Mass2_3)##
## Call:
## pgls(formula = WinFreq ~ Mass2_3, data = video_comparative, lambda = "ML")
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.031340 -0.010905 -0.002889 0.022627 0.064360
##
## Branch length transformations:
##
## kappa [Fix] : 1.000
## lambda [ ML] : 0.558
## lower bound : 0.000, p = 0.21225
## upper bound : 1.000, p = 0.16822
## 95.0% CI : (NA, NA)
## delta [Fix] : 1.000
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.2890066 0.1631528 1.7714 0.09825 .
## Mass2_3 -0.0025029 0.0059583 -0.4201 0.68081
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.02923 on 14 degrees of freedom
## Multiple R-squared: 0.01245, Adjusted R-squared: -0.05809
## F-statistic: 0.1765 on 1 and 14 DF, p-value: 0.6808# Check the distribution of the residuals:
par(mfrow=c(2,2))
plot(live_pgls_CBI_Mass)
par(mfrow=c(2,2))
plot(live_pgls_CBI_Mass2_3)
par(mfrow=c(2,2))
plot(live_pgls_WinFreq_Mass)
par(mfrow=c(2,2))
plot(live_pgls_WinFreq_Mass2_3)
par(mfrow=c(2,2))
plot(live_pgls_CBI_Mass_nosp)
par(mfrow=c(2,2))
plot(live_pgls_CBI_Mass2_3_nosp)
par(mfrow=c(2,2))
plot(live_pgls_WinFreq_Mass_nosp)
par(mfrow=c(2,2))
plot(live_pgls_WinFreq_Mass2_3_nosp)
par(mfrow=c(2,2))
plot(video_pgls_CBI_Mass)
par(mfrow=c(2,2))
plot(video_pgls_CBI_Mass2_3)
par(mfrow=c(2,2))
plot(video_pgls_WinFreq_Mass)
par(mfrow=c(2,2))
plot(video_pgls_WinFreq_Mass2_3)
Generate dataset-specific versions of Figure 1:
png("/Users/David/Documents/College/2017–18/Winter 2018/FBQ paper/Behav Ecol revisions/Supp_Fig1_live.png", width = 3600, height = 3600, pointsize = 96)
layout(matrix(c(seq(1,4,1)), nrow=2, ncol=2, byrow=TRUE))
par(mar=c(5, 4, 2, 2))
plot(sheet_phylo_test$Mass, sheet_phylo_test$CBI, xlab = "Body mass (g)", ylab = "CBI", main = "Isometric scaling", pch = 16, col = ifelse(sheet_phylo_test$Mass == sheet_phylo_test$Mass[sheet_phylo_test$Species == "Spotted Porcupinefish"], "gray70", "black"))
title("(a)", line = -1, adj = 0, outer = TRUE)
title("(b)", line = -1, outer = TRUE)
title("(c)", line = -20, adj = 0, outer = TRUE)
title("(d)", line = -20, outer = TRUE)
abline(live_pgls_CBI_Mass, lwd = 5)
abline(live_pgls_CBI_Mass_nosp, lwd = 5, col = "gray70")
lambda <- live_pgls_CBI_Mass$param[2]
pseudor2 <- summary(live_pgls_CBI_Mass)$r.squared
pvalue <- summary(live_pgls_CBI_Mass)$coef[2,4]
lambda_nosp <- live_pgls_CBI_Mass_nosp$param[2]
pseudor2_nosp <- summary(live_pgls_CBI_Mass_nosp)$r.squared
pvalue_nosp <- summary(live_pgls_CBI_Mass_nosp)$coef[2,4]
text(600, 2.95, cex = 0.8, bquote(lambda == .(spec_dec(lambda, 3))), pos = 4)
text(600, 2.75, cex = 0.8, bquote(italic(R)^2 == .(spec_dec(pseudor2, 3))), pos = 4)
text(600, 2.5, cex = 0.8, bquote(italic(p) == .(spec_dec(pvalue, 3))), pos = 4)
text(600, 1.6, cex = 0.8, bquote(lambda == .(spec_dec(lambda_nosp, 3))), col = "gray70", pos = 4)
text(600, 1.4, cex = 0.8, bquote(italic(R)^2 == .(spec_dec(pseudor2_nosp, 3))), col = "gray70", pos = 4)
text(600, 1.15, cex = 0.8, bquote(italic(p) == .(spec_dec(pvalue_nosp, 3))), col = "gray70", pos = 4)
par(mar=c(5, 4, 2, 2))
plot(sheet_phylo_test$Mass2_3, sheet_phylo_test$CBI, xlab = "Body mass raised to 2/3", ylab = "", main = "Allometric scaling", pch = 16, col = ifelse(sheet_phylo_test$Mass == sheet_phylo_test$Mass[sheet_phylo_test$Species == "Spotted Porcupinefish"], "gray70", "black"))
abline(live_pgls_CBI_Mass2_3, lwd = 5)
abline(live_pgls_CBI_Mass2_3_nosp, lwd = 5, col = "gray70")
lambda <- live_pgls_CBI_Mass2_3$param[2]
pseudor2 <- summary(live_pgls_CBI_Mass2_3)$r.squared
pvalue <- summary(live_pgls_CBI_Mass2_3)$coef[2,4]
lambda_nosp <- live_pgls_CBI_Mass2_3_nosp$param[2]
pseudor2_nosp <- summary(live_pgls_CBI_Mass2_3_nosp)$r.squared
pvalue_nosp <- summary(live_pgls_CBI_Mass2_3_nosp)$coef[2,4]
text(58, 2.75, cex = 0.8, bquote(lambda == .(spec_dec(lambda, 3))), pos = 4)
text(58, 2.55, cex = 0.8, bquote(italic(R)^2 == .(spec_dec(pseudor2, 3))), pos = 4)
text(58, 2.3, cex = 0.8, bquote(italic(p) == .(spec_dec(pvalue, 3))), pos = 4)
text(58, 0.6, cex = 0.8, bquote(lambda == .(spec_dec(lambda_nosp, 3))), col = "gray70", pos = 4)
text(58, 0.4, cex = 0.8, bquote(italic(R)^2 == .(spec_dec(pseudor2_nosp, 3))), col = "gray70", pos = 4)
text(58, 0.15, cex = 0.8, bquote(italic(p) == .(spec_dec(pvalue_nosp, 3))), col = "gray70", pos = 4)
par(mar=c(5, 4, 2, 2))
plot(sheet_phylo_test$Mass, sheet_phylo_test$WinFreq, xlab = "Body mass (g)", ylab = "Frequency of winning", pch = 16, col = ifelse(sheet_phylo_test$Mass == sheet_phylo_test$Mass[sheet_phylo_test$Species == "Spotted Porcupinefish"], "gray70", "black"))
abline(live_pgls_WinFreq_Mass, lwd = 5)
abline(live_pgls_WinFreq_Mass_nosp, lwd = 5, col = "gray70")
lambda <- live_pgls_WinFreq_Mass$param[2]
pseudor2 <- summary(live_pgls_WinFreq_Mass)$r.squared
pvalue <- summary(live_pgls_WinFreq_Mass)$coef[2,4]
lambda_nosp <- live_pgls_WinFreq_Mass_nosp$param[2]
pseudor2_nosp <- summary(live_pgls_WinFreq_Mass_nosp)$r.squared
pvalue_nosp <- summary(live_pgls_WinFreq_Mass_nosp)$coef[2,4]
text(600, 0.9, cex = 0.8, bquote(lambda == .(spec_dec(lambda, 3))), pos = 4)
text(600, 0.9 - 1.2/18, cex = 0.8, bquote(italic(R)^2 == .(spec_dec(pseudor2, 3))), pos = 4)
text(600, 0.9 - 1.2/8, cex = 0.8, bquote(italic(p) == .(spec_dec(pvalue, 3))), pos = 4)
text(600, 0.35, cex = 0.8, bquote(lambda == .(spec_dec(lambda_nosp, 3))), col = "gray70", pos = 4)
text(600, 0.35 - 1.2/18, cex = 0.8, bquote(italic(R)^2 == .(spec_dec(pseudor2_nosp, 3))), col = "gray70", pos = 4)
text(600, 0.35 - 1.2/8, cex = 0.8, bquote(italic(p) == .(spec_dec(pvalue_nosp, 3))), col = "gray70", pos = 4)
par(mar=c(5, 4, 2, 2))
plot(sheet_phylo_test$Mass2_3, sheet_phylo_test$WinFreq, xlab = "Body mass raised to 2/3", ylab = "", pch = 16, col = ifelse(sheet_phylo_test$Mass == sheet_phylo_test$Mass[sheet_phylo_test$Species == "Spotted Porcupinefish"], "gray70", "black"))
abline(live_pgls_WinFreq_Mass2_3, lwd = 5)
abline(live_pgls_WinFreq_Mass2_3_nosp, lwd = 5, col = "gray70")
lambda <- live_pgls_WinFreq_Mass2_3$param[2]
pseudor2 <- summary(live_pgls_WinFreq_Mass2_3)$r.squared
pvalue <- summary(live_pgls_WinFreq_Mass2_3)$coef[2,4]
lambda_nosp <- live_pgls_WinFreq_Mass2_3_nosp$param[2]
pseudor2_nosp <- summary(live_pgls_WinFreq_Mass2_3_nosp)$r.squared
pvalue_nosp <- summary(live_pgls_WinFreq_Mass2_3_nosp)$coef[2,4]
text(58, 0.5, cex = 0.8, bquote(lambda == .(spec_dec(lambda, 3))), pos = 4)
text(58, 0.5 - 1.2/18, cex = 0.8, bquote(italic(R)^2 == .(spec_dec(pseudor2, 3))), pos = 4)
text(58, 0.5 - 1.2/8, cex = 0.8, bquote(italic(p) == .(spec_dec(pvalue, 3))), pos = 4)
text(58, 0.2, cex = 0.8, bquote(lambda == .(spec_dec(lambda_nosp, 3))), col = "gray70", pos = 4)
text(58, 0.2 - 1.2/18, cex = 0.8, bquote(italic(R)^2 == .(spec_dec(pseudor2_nosp, 3))), col = "gray70", pos = 4)
text(58, 0.2 - 1.2/8, cex = 0.8, bquote(italic(p) == .(spec_dec(pvalue_nosp, 3))), col = "gray70", pos = 4)
dev.off()
png("/Users/David/Documents/College/2017–18/Winter 2018/FBQ paper/Behav Ecol revisions/Supp_Fig1_video.png", width = 3600, height = 3600, pointsize = 96)
layout(matrix(c(seq(1,4,1)), nrow=2, ncol=2, byrow=TRUE))
par(mar=c(5, 4, 2, 2))
plot(scored_phylo_test$Mass, scored_phylo_test$CBI, xlab = "Body mass (g)", ylab = "CBI", main = "Isometric scaling", pch = 16, col = "black")
title("(a)", line = -1, adj = 0, outer = TRUE)
title("(b)", line = -1, outer = TRUE)
title("(c)", line = -20, adj = 0, outer = TRUE)
title("(d)", line = -20, outer = TRUE)
abline(video_pgls_CBI_Mass, lwd = 5)
lambda <- video_pgls_CBI_Mass$param[2]
pseudor2 <- summary(video_pgls_CBI_Mass)$r.squared
pvalue <- summary(video_pgls_CBI_Mass)$coef[2,4]
text(150, 2.5, cex = 0.8, bquote(lambda == .(spec_dec(lambda, 3))), pos = 4)
text(150, 2.3, cex = 0.8, bquote(italic(R)^2 == .(spec_dec(pseudor2, 3))), pos = 4)
text(150, 2.05, cex = 0.8, bquote(italic(p) == .(spec_dec(pvalue, 3))), pos = 4)
par(mar=c(5, 4, 2, 2))
plot(scored_phylo_test$Mass2_3, scored_phylo_test$CBI, xlab = "Body mass raised to 2/3", ylab = "", main = "Allometric scaling", pch = 16, col = "black")
abline(video_pgls_CBI_Mass2_3, lwd = 5)
lambda <- video_pgls_CBI_Mass2_3$param[2]
pseudor2 <- summary(video_pgls_CBI_Mass2_3)$r.squared
pvalue <- summary(video_pgls_CBI_Mass2_3)$coef[2,4]
text(25, 2.5, cex = 0.8, bquote(lambda == .(spec_dec(lambda, 3))), pos = 4)
text(25, 2.3, cex = 0.8, bquote(italic(R)^2 == .(spec_dec(pseudor2, 3))), pos = 4)
text(25, 2.05, cex = 0.8, bquote(italic(p) == .(spec_dec(pvalue, 3))), pos = 4)
par(mar=c(5, 4, 2, 2))
plot(scored_phylo_test$Mass, scored_phylo_test$WinFreq, xlab = "Body mass (g)", ylab = "Frequency of winning", pch = 16, col = "black")
abline(video_pgls_WinFreq_Mass, lwd = 5)
lambda <- video_pgls_WinFreq_Mass$param[2]
pseudor2 <- summary(video_pgls_WinFreq_Mass)$r.squared
pvalue <- summary(video_pgls_WinFreq_Mass)$coef[2,4]
text(150, 0.55, cex = 0.8, bquote(lambda == .(spec_dec(lambda, 3))), pos = 4)
text(150, 0.55 - 1.2/18, cex = 0.8, bquote(italic(R)^2 == .(spec_dec(pseudor2, 3))), pos = 4)
text(150, 0.55 - 1.2/8, cex = 0.8, bquote(italic(p) == .(spec_dec(pvalue, 3))), pos = 4)
par(mar=c(5, 4, 2, 2))
plot(scored_phylo_test$Mass2_3, scored_phylo_test$WinFreq, xlab = "Body mass raised to 2/3", ylab = "", pch = 16, col = "black")
abline(video_pgls_WinFreq_Mass2_3, lwd = 5)
lambda <- video_pgls_WinFreq_Mass2_3$param[2]
pseudor2 <- summary(video_pgls_WinFreq_Mass2_3)$r.squared
pvalue <- summary(video_pgls_WinFreq_Mass2_3)$coef[2,4]
text(25, 0.55, cex = 0.8, bquote(lambda == .(spec_dec(lambda, 3))), pos = 4)
text(25, 0.55 - 1.2/18, cex = 0.8, bquote(italic(R)^2 == .(spec_dec(pseudor2, 3))), pos = 4)
text(25, 0.55 - 1.2/8, cex = 0.8, bquote(italic(p) == .(spec_dec(pvalue, 3))), pos = 4)
dev.off()Fit Lanchester’s laws as well as the body size-only and group size-only models to the live-only and video-only datasets:
live_group <- sheet4[sheet4$Group.A != sheet4$Group.B,]
video_group <- scored4[scored4$Group.A != scored4$Group.B,]
# Exclude cornetfish:
live_group <- live_group[live_group$Species.A != "Cornetfish" &
live_group$Species.B != "Cornetfish",]
# Number of data points:
nrow(live_group)## [1] 59nrow(video_group)## [1] 24# Get body size info and use it to test the four models:
lanchester_on_frames <- function(dataframe) {
n <- ncol(dataframe)
dataframe$Group.A <- as.numeric(dataframe$Group.A)
dataframe$Group.B <- as.numeric(dataframe$Group.B)
for(i in 1:nrow(dataframe)) {
dataframe[i, n+1] <- lengthmass$Mass[lengthmass$Species == as.character(dataframe$Species.A[i])]
dataframe[i, n+2] <- lengthmass$Mass[lengthmass$Species == as.character(dataframe$Species.B[i])]
dataframe[i, n+3] <- lengthmass$Mass2_3[lengthmass$Species == as.character(dataframe$Species.A[i])]
dataframe[i, n+4] <- lengthmass$Mass2_3[lengthmass$Species == as.character(dataframe$Species.B[i])]
if (dataframe$Group.A[i]*dataframe[i, n+1] > dataframe$Group.B[i]*dataframe[i, n+2]) {
dataframe[i, n+5] <- dataframe$Species.A[i]
} else {
dataframe[i, n+5] <- dataframe$Species.B[i]
}
if ((dataframe$Group.A[i]^2)*dataframe[i, n+1] > (dataframe$Group.B[i]^2)*dataframe[i, n+2]) {
dataframe[i, n+6] <- dataframe$Species.A[i]
} else {
dataframe[i, n+6] <- dataframe$Species.B[i]
}
if (identical(as.character(dataframe$Winner[i]), as.character(dataframe[i, n+5])) == TRUE) {
dataframe[i, n+7] <- 1
} else {
dataframe[i, n+7] <- 0
}
if (identical(as.character(dataframe$Winner[i]), as.character(dataframe[i, n+6])) == TRUE) {
dataframe[i, n+8] <- 1
} else {
dataframe[i, n+8] <- 0
}
if (dataframe$Group.A[i]*dataframe[i, n+3] > dataframe$Group.B[i]*dataframe[i, n+4]) {
dataframe[i, n+9] <- dataframe$Species.A[i]
} else {
dataframe[i, n+9] <- dataframe$Species.B[i]
}
if ((dataframe$Group.A[i]^2)*dataframe[i, n+3] > (dataframe$Group.B[i]^2)*dataframe[i, n+4]) {
dataframe[i, n+10] <- dataframe$Species.A[i]
} else {
dataframe[i, n+10] <- dataframe$Species.B[i]
}
if (identical(as.character(dataframe$Winner[i]), as.character(dataframe[i, n+9])) == TRUE) {
dataframe[i, n+11] <- 1
} else {
dataframe[i, n+11] <- 0
}
if (identical(as.character(dataframe$Winner[i]), as.character(dataframe[i, n+10])) == TRUE) {
dataframe[i, n+12] <- 1
} else {
dataframe[i, n+12] <- 0
}
if (dataframe$Group.A[i] > dataframe$Group.B[i]) {
dataframe[i, n+13] <- dataframe$Species.A[i]
} else {
dataframe[i, n+13] <- dataframe$Species.B[i]
}
if (dataframe[i, n+1] > dataframe[i, n+2]) {
dataframe[i, n+14] <- dataframe$Species.A[i]
} else {
dataframe[i, n+14] <- dataframe$Species.B[i]
}
if (identical(as.character(dataframe$Winner[i]), as.character(dataframe[i, n+13])) == TRUE) {
dataframe[i, n+15] <- 1
} else {
dataframe[i, n+15] <- 0
}
if (identical(as.character(dataframe$Winner[i]), as.character(dataframe[i, n+14])) == TRUE) {
dataframe[i, n+16] <- 1
} else {
dataframe[i, n+16] <- 0
}
}
colnames(dataframe)[(n+1):(n+16)] <- c("BodySize.A", "BodySize.B", "BodySize2_3.A",
"BodySize2_3.B", "Linear.Law.Prediction",
"Square.Law.Prediction", "Linear.True",
"Square.True", "Linear.Law2_3.Prediction",
"Square.Law2_3.Prediction", "Linear2_3.True",
"Square2_3.True", "Group.Size.Only.Prediction",
"Body.Size.Only.Prediction",
"Group.Size.Only.True", "Body.Size.Only.True")
return(dataframe)
}
live_group <- lanchester_on_frames(live_group)
video_group <- lanchester_on_frames(video_group)
# Proportion of successful predictions:
sum(live_group$Linear.True)/nrow(live_group)## [1] 0.6949153sum(live_group$Square.True)/nrow(live_group)## [1] 0.5254237sum(live_group$Linear2_3.True)/nrow(live_group)## [1] 0.6101695sum(live_group$Square2_3.True)/nrow(live_group)## [1] 0.4237288sum(live_group$Group.Size.Only.True)/nrow(live_group)## [1] 0.2542373sum(live_group$Body.Size.Only.True)/nrow(live_group)## [1] 0.6779661sum(video_group$Linear.True)/nrow(video_group)## [1] 0.375sum(video_group$Square.True)/nrow(video_group)## [1] 0.08333333sum(video_group$Linear2_3.True)/nrow(video_group)## [1] 0.1666667sum(video_group$Square2_3.True)/nrow(video_group)## [1] 0.04166667sum(video_group$Group.Size.Only.True)/nrow(video_group)## [1] 0sum(video_group$Body.Size.Only.True)/nrow(video_group)## [1] 0.75# Test significance:
dbinom(sum(live_group$Linear.True), size = nrow(live_group), prob = 1/2) ## [1] 0.001123269dbinom(sum(live_group$Square.True), size = nrow(live_group), prob = 1/2)## [1] 0.09596023dbinom(sum(live_group$Linear2_3.True), size = nrow(live_group), prob = 1/2) ## [1] 0.02501637dbinom(sum(live_group$Square2_3.True), size = nrow(live_group), prob = 1/2)## [1] 0.05253438dbinom(sum(live_group$Group.Size.Only.True), size = nrow(live_group), prob = 1/2) ## [1] 6.920778e-05dbinom(sum(live_group$Body.Size.Only.True), size = nrow(live_group), prob = 1/2)## [1] 0.002423897dbinom(sum(video_group$Linear.True), size = nrow(video_group), prob = 1/2) ## [1] 0.07793331dbinom(sum(video_group$Square.True), size = nrow(video_group), prob = 1/2)## [1] 1.645088e-05dbinom(sum(video_group$Linear2_3.True), size = nrow(video_group), prob = 1/2) ## [1] 0.000633359dbinom(sum(video_group$Square2_3.True), size = nrow(video_group), prob = 1/2)## [1] 1.430511e-06dbinom(sum(video_group$Group.Size.Only.True), size = nrow(video_group), prob = 1/2) ## [1] 5.960464e-08dbinom(sum(video_group$Body.Size.Only.True), size = nrow(video_group), prob = 1/2)## [1] 0.008022547How many interactions involving damselfish do we have in the live-only dataset and how many in the video-only dataset?
nrow(live_group[live_group$Species.A == "Damselfish" | live_group$Species.B == "Damselfish",])/nrow(live_group)## [1] 0.3559322nrow(video_group[video_group$Species.A == "Damselfish" | video_group$Species.B == "Damselfish",])/nrow(video_group)## [1] 0.7916667References
- Clutton-Brock TH, Albon SD, Gibson RM, Guinness FE 1979 The logical stag: adaptive aspects of fighting in red deer (Cervus elaphus L.). Anim Behav 27: 211–225
- Clutton-Brock TH, Guinness FE, Albon SD 1982 Red Deer: Behaviour and Ecology of Two Sexes. Chicago, IL: Univ of Chicago Press
- Randall JE 2005 Reef and Shore Fishes of the South Pacific: New Caledonia to Tahiti and the Pitcairn Islands. Honolulu, HI: Univ of Hawaii Press
- Schmid VS, de Vries H. 2013. Finding a dominance order most consistent with a linear hierarchy: an improved algorithm for the I&SI method. Anim Behav 86: 1097–1105
- Shelley EL, Tanaka MYU, Ratnathicam AR, Blumstein DT. 2004. Can Lanchester’s laws help explain interspecific dominance in birds? Condor 106: 395–400